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BISTATIC  SCATTERER  AND  ANTENNA  IMAGING 

Steve  Inge  (in-j) 

Broadcast  Communications  Ltd. 

P.O.  Box  98,  Wellington,  New  Zealand 


ABSTRACT:  Two  computing  techniques  to  create 
an  image  of  the  radiating  centres  on  an  antenna  or 
scatterer  using  Fourier  optics  is  presented. 


Steve  Inge  died  on  October  24,  2002  at  the  age  of  61  after  this 
article  was  submitted  for  publication.  The  reviewers  accepted 
this  article  subject  to  revisions.  Since  that  is  not  possible,  this 
article  is  published  as  submitted. 

The  late  Ian  McEnnis,  Chief  Antenna  Engineer  at  Broadcast 
Communications  Limited  in  New  Zealand,  and  one  of  Steve*s 
colleagues,  notified  us  of  his  death  and  had  the  following 
remembrances  of  Steve,  which  we  print  below. 

Steve  had  been  with  our  organization  for  42  years.  His  immense 
intellect  enabled  him  to  become  a  **gurtt*’  in  any  area  he  focused 
on.  Most  of  BCL’s  in  house  engineering  sofhvare  has  been 
developed  by  Steve.  BCLIPPS  (BCL*s  Interactive  Planning  and 
Propagation  Software)  was  Steve’s  crowning.  It  has  formed  the 
basis  of  our  coverage  work  for  the  past  10  years  and  has  enabled 
us  to  be  a  leader  in  this  area.  Steve  was  the  engineer  who  solved 
the  most  difficult  of  problems  and  he  was  the  one  his  colleagues 
turned  to  when  they  were  having  difficulty.  His  mathematical 
ability  was  awesome  and  once  he  got  his  teeth  into  a  problem  he 
would  not  let  it  go  until  he  had  it  sorted.  There  is  nothing  Steve 
liked  more  than  to  share  his  knowledge  with  others.  Having  said 
all  that,  in  the  end  Steve  was  Just  a  ’’good  bloke**  who  will  be 
sorely  missed  by  many. 


Introduction 

Bistatic  k-space  (BSKS)  imaging  from  antenna 
currents  was  introduced  by  John  Shaeffer  et  al  (Ref 
1).  It  is  a  technique  to  graphically  show  the  active 


radiating  centers  on  a  conducting  body.  The  currents 
on  the  body  are  found  using  either  a  computer 
program  or  measurements,  for  a  particular  excitation. 
The  excitation  may  be  a  plane  wave  being  scattered 
by  the  body,  or  sources  on  the  body  for  an  antenna. 
This  paper  describes  two  simple  approaches  to 
imaging  a  radiator. 

Brief  Overview 

The  radiation  emanating  from  any  antenna  (or 
scatterer)  may  be  visualised  by  simple  operations 
upon  the  antenna  currents.  First  choose  the  far  field 
direction  from  which  to  look  at  the  antenna.  Create  a 
plane  surface  at  the  antenna,  normal  to  the  look 
direction.  Project  all  of  the  antenna  currents  onto  this 
surface,  adjusting  the  phase  as  necessary,  to  maintain 
the  same  phases  and  polarisation  in  the  far  field. 
Filter  the  current  image  on  the  surface  by  removing 
all  spatial  frequencies  above  one  cycle  per 
wavelength.  Display  the  current  intensiQr  as 
brightness  in  a  2D  graphic  image, 'which  will  show 
which  parts  of  the  antenna  are  "glowing"  when 
viewed  from  this  direction.  This  is  repeated  for  both 
polarisations,  which  may  be  combined  into  a  single 
total  radiation  image  if  desired. 

Imaging  theory 

How  does  it  work?  By  flattening  the  currents  on  to  a 
surface  we  create  a  two  dimensional  current  plane.  If 
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we  now  take  the  Fourier  transform  (FT)  of  this 
current  surface,  we  create  a  new  description  of  the 
currents.  The  FT  creates  a  set  of  sinusoidal  waves 
equivalent  to  the  original  currents.  Since  the  currents 
occupied  a  finite  area  on  the  surface,  the  FT  range  is 
infinite.  Each  point  on  the  FT  surface  represents  a 
particular  uniform  linear  phase  advance,  and  thus 
represents  a  plane  wave  propagating  in  a  single 
direction.  This  is  known  as  a  plane  wave  spectrum  of 
the  radiation.  The  Fx,  Fy  positions  on  the  spectrum 
are  the  direction  cosines  of  the  plane  wave  direction. 
The  z  component  of  the  radiation  is 

Fz  =  yj\-(Fx^  +  Fy^)  The  centre  of  the  plane  (0,0)  is 

constant  phase,  which  creates  a  wavefront 
propagating  toward  the  look  direction.  When 
Fx^  +  Fy^  the  wave  direction  is  normal  to  the 
look  direction.  All  points  on  this  unit  circle 
correspond  to  a  phase  advance  of  one  cycle  per 
wavelength.  The  spectrum  outside  this  circle  gives 
direction  cosines  >  I,  indicating  imaginary  angles. 
This  energy  does  not  radiate,  but  represents 
evanescent  waves  which  eventually  return  to  the 
antenna  as  reactive  energy.  (Ref  3) 

The  antenna  image  is  created  by  projecting  the 
currents  onto  a  normal  plane,  and  removing  spatial 
frequencies  greater  than  one  cycle  per  wavelength. 
The  latter  process  may  be  achieved  by  setting  the  FT 
spectrum  to  zero  at  all  points  outside  the  one  cycle 
per  wavelength  circle  and  inverse  FT,  or  equivalently 
by  convolving  the  current  surface  with  the 
J^{np)/{Kp)  point  spread  function  (PSF),  where 

7]  0  is  the  Bessel  function  of  the  first  kind  and  order 
one,  and  rho  is  the  radius.  This  point  spread  function 
is  known  from  optics  as  the  Airy  disk. 


The  image  is  a  view  of  the  far  field  antenna  radiation 
toward  the  look  direction.  For  scatterers,  the  look 
direction  may  be  the  same  as  the  incident 
illumination,  which  is  the  monostatic  (single  antenna) 
radar  case.  Similarly,  when  the  look  direction  is 
different  to  the  illumination  direction,  we  have  the 
bistatic  (separate  antennas)  radar  case. 

Imaging  process 

1.  Given  the  locations  and  currents  of  wire  segments, 
from  a  NEC  (Ref  4)  analysis  for  example,  calculate 
the  contribution  to  far  field  (Efar)  for  each 
segment,  in  a  chosen  (theta,  phi)  radiation 
direction.  The  two  polarisations  are  treated 
separately.  The  contribution  from  a  single  segment 
is 

A£^=y30A/expOAR.S){(R«AL)R-AL}*p  (1) 
where  p  is  the  chosen  polarisation  unit  vector,  I  is 
the  current,  AL  is  the  segment  length  vector,  R  is 

the  far  field  direction  unit  vector,  S  is  the 
coordinate  vector  of  the  segment  centre  and  k  is  the 
wave  number  Itc/Z.  Note  that  the  1/r  distance 
factor  is  omitted,  making  the  units  in  volts.  This  is 
the  same  operation  as  calculating  the  far  field 
radiation  pattern,  without  summing  the  voltages." 
The  code  can  be  extracted  from  a  suitable  antenna 
program  far  field  radiation  pattern  routine. 

2.  Each  segment  is  replaced  with  an  isotropic  source 
in  the  same  location,  with  the  magnitude  and  phase 
of  the  Efar  contribution  to  the  chosen  polarisation. 
By  placing  this  light  source  at  each  point,  we  create 
an  image,  viewable  from  any  direction.  At  this 
point  all  the  detail  can  be  seen,  so  that  wires  with 
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rapidly  varying  phase  from  a  particular  viewpoint 
are  visible.  This  does  not  agree  with  radiation 
reality. 

3.  Filter  the  image  to  remove  spatial  frequencies 
above  one  cycle  per  wavelength.  This  physical 
limit  is  due  solely  to  the  wave  nature  of  light.  Note 
that  even  the  tiniest  point  source  appears  as  a 
circular  blob,  surrounded  by  faint  rings,  known  as 
the  Airy  disk.  The  Aiiy  disk  has  the  form  of 
Jj(M)/(t/),  where  is  the  Bessel  function  of 
the  first  kind  and  order  1.  The  resolution  limit  is 
the  Raleigh  criterion  of  0.61  wavelength,  which  is 
where  the  maximum  of  one  Airy  disk  overlaps  with 
the  first  minimum  of  another. 

Starting  with  the  Efar  elements,  we  Fast  Fourier 
Transform  (FFT)  the  (ID,  2D,  or  3D)  space, 
remove  all  spatial  frequencies  above  one  cycle  per 
wavelength,  and  inverse  FFT.  The  wanted  spatial 
frequencies  are  found  in  the  line,  circle  or  sphere 
sectors  in  the  comers  of  the  transformed 
matrix,  the  remainder  are  set  to  zero. 

4.  Now  the  image  appears  as  it  would  in  our 
“microwave  microscope”.  Note  that  optical 
microscopes  have  the  same  limit  to  magnification, 
where  a  point  is  reached  where  more  magnificati<m^ 
just  gives  more  blur.  Hence  the  move  to  higher 
frequencies  and  electron  microscopes.  These 
images  may  be  named  radio  wave  photographs, 
since  they  show  the  luminosity  of  each  part  of  the 
antenna  or  scatterer. 

Example  Consider  a  long  wire  with  a  sinusoidal 

standing  or  traveling  wave  current  pattern.  If  we  view 


the  currents  from  broadside,  we  can  see  the 
alternating  pattern.  However  if  we  apply  the  blur 
filter,  we  find  that  at  all  internal  points  along  the 
wire,  the  cyclic  currents  integrate  to  zero.  Only  at  the 
ends  is  the  last  current  pattern  not  cancelled  by  its 
neighbour,  so  we  see  a  bright  spot  centred  on  each 
wire  end.  This  illustrates  the  maxim  "radiation  is  due 
to  the  acceleration  of  charge".  See  the  first  example 
below. 

Taper  windows 

Tapering  of  the  spatial  window  illumination  reduces 
both  resolution  and  sidelobe  levels.  The  taper 

function  used  is  where  rho  is  the  fraction 

of  the  window  radius,  and  p  is  the  taper  exponent.  In 
2D,  the  resulting  blur  function  form  is 

(u)  / giving  both  reduced  resolution  and 

sidelobe  levels  as  p  increases.  />  =  0  for  uniform 
illumination,  p  =  1  for  parabolic  taper,  etc. 

i 

Example  radio  wave  photographs 

The  grey  scale  is  in  decibels.  The  number  at  top  left 
of  each  image  is  the  maximum  "brightness"  across 
the  image.  The  look  direction  for  all  the  images  is 
from  the  zenith,  ^  =  0,  (p  =1®.  The  small-fp  offset  is 
to  ensure  that  the  polarisation  is  properly  defined.  All 
these  images  used  NEC-2  to  calculate  the  current 
distributions.  Unless  specified,  horizontal  (Ephi) 
polarisation  is  used.  The  faint  grid  lines  are  spaced 
one  wavelength  apart.  A  multidimensional  mixed 
radix  FFT  is  used  (Ref  5),  with  each  image  as 
200x200  pixels.  In  the  largest  image  spaces  the  pixel 
size  is  )J20. 


Tab!e  1.  Window  Taper  Performance 
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Radio  Camera 

An  alternative  imaging  process  is  to  emulate  a 
camera.  In  this  process  a  plane  is  created  in  the  far 
field  of  the  antenna/scatterer,  and  the  E  field  is 
calculated  at  many  points  on  the  plane  (usually  in  a 
rectangular  grid).  This  array  of  complex  field  values 
is  a  numerical  hologram.  A  reference  wave  is  not 
needed,  as  the  phase  is  recorded  in  the  complex  E 
field  value.  This  numerical  hologram  may  be 
reconstructed  using  a  camera  lens,  or  equivalently  by 
taking  the  FT  of  the  2D  complex  field  values.  The 
image  is  formed  at  one  focal  length  behind  the 
camera  lens.  The  radiation  integral  has  already 
removed  the  evanescent  energy,  so  it  is  not  necessary 
to  filter  the  image  any  further.  Note  that  the  point 
spread  function  may  be  larger  than  the  theoretical 
Airy  disk,  since  the  finite  lens  aperture  (or  grid  size) 
can  only  worsen  the  resolution.  In  microscopy  this 
reduction  is  known  as  the  “numerical  aperture”. 
Using  NEC  for  example,  we  can  evaluate  the  far  field 
over  a  rectangular  grid  of  points  on  a  plane  normal  to 
the  look  angle.  The  difference  from  the  previous 
imaging  process  is  that  the  plane  is  in  the  far  field. 
The  antenna/scatterer  should  be  rotated  (with  a  GM 
card)  so  that  the  desired  look  angle  is  moved  to  the  x 
axis.  The  near  field  command  (NE)  is  used  to  create  a 
planar  grid  of  field  values  in  the  YZ  plane.  There  are 
still  two -polarisations  required.  After  gathering  the. 
output  values  into  a  suitable  matrix  form,  the  FT 
image  may  be  formed  using  Matlab  or  Mathematica 
etc. 
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Conclusion 

Two  practical  radiator  imaging  techniques  have  been 
described.  These  are  thought  to  be  new,  but 
similarities  to  near  field  measuring  techniques  and 
holography  suggests  that  they  have  simply  been 
rediscovered.  However,  it  is  timely  to  bring  imaging 
to  the  attention  of  CEMists  as  a  useful  tool  in 
understanding  radiation.  The  author  has  used  the  first 
technique  to  examine  currents  induced  in  a  support 
tower,  which  caused  an  unwanted  null  in  the 
radiation  pattern.  Viewing  the  currents  only  confused 
the  situation,  since  most  of  the  currents  did  not 
radiate  in  the  null  direction. 
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Abstract:  Comprehensive  shielding  in  modern  electronic 
equipment  may  lead  to  resonant  behaviour  within  the 
equipment  enclosure.  This  paper  presents  a  method  for 
optimising  the  placement  of  sources  of  electromagnetic 
(EM)  energy  and  susceptors  to  this  EM  energy  within  an 
enclosed  resonant  cavity.  The  source  and  susceptor  2tre 
placed  on  a  dielectric  slab  within  a  perfectly  conducting 
enclosure  to  reduce  the  level  of  EM  coupling  between  the 
two.  Optimisation  is  facilitated  using  a  genetic  algorithm 
coupled  with  a  finite-difference  time-domain  solver.  Re¬ 
sults  are  presented  for  single  objective  optimisation  and 
multi-objective  optimisation  cases. 

1  Introduction 

The  explosive  increase  in  the  use  of  electronic  equipment 
in  the  information  age  has  led  to  electromagnetic  com¬ 
patibility  (EMC)  and  electromagnetic  interference  (EMI) 
issues  becoming  more  important  to  designers  of  electronic 
equipment.  These  issues  must  be  considered  at  the  time 
of  design  and  not  once  designs  are  at  the  prototype  stage. 
Higher  clock  speeds  and  faster  switching  transitions  lead 
to  greater  levels  of  electromagnetic  emissions,  while  higher 
component  integration  and  lower  power  demands  lead  to 
greater  sensitivity.  So  at  once  systems  are  becoming  more 
prone  to  generating  and  also  more  sensitive  towards  EMI. 

An  area  of  interest  to  the  authors  is  electronic  enclosures 
and  the  placement  of  devices  inside  these  enclosures.  Mod¬ 
ern  electronic  items  have  many  different  sources  of  EMI. 
These  sources  are  often  placed  inside  a  rectangular  shaped 
metal  box  to  limit  the  EM  emissions  from  the  equipment. 
Care  is  taken  to  ensure  the  shielding  is  as  comprehensive  as 
possible  thus  apertures  are  kept  to  a  minimum  and  covered 
in  metallic  blanking  plates.  Such  a  shielded,  rectangular 
enclosure  has  a  good  chance  of  forming  a  resonant  cavity 
so  that  field  strengths  generated  by  the  circuit  may  well 
be  enhanced  once  the  circuit  is  mounted  within  the  casing. 
Normally  the  circuit  may  have  a  number  of  sources  within 
the  enclosure,  all  of  which  add  to  the  EM  fields  within; 
changing  the  positions  of  the  various  sources  will  change 
the  amount  of  resonance  and  constructive  interference,  [1]. 
There  are  also  likely  to  be  a  number  of  susceptors  at  mul¬ 
tiple  locations  throughout  the  circuit. 


There  is  obviously  an  optimal  size  and  shape  of  enclo¬ 
sure  and  component  layout  however,  exhaustively  placing 
components  inside  enclosures  and  then  computing  the  re¬ 
sultant  fields  is  a  task  that  that  would  require  a  massive 
undertaking  on  behalf  of  the  designer.  A  far  better  ap¬ 
proach  to  component  placement  is  to  use  some  kind  of  op¬ 
timisation  method  to  place  the  components  to  achieve  a 
desired  radiation  level.  This  paper  describes  the  novel  use 
of  a  genetic  algorithm  (GA)  and  a  finite-difference  time- 
domain  (FDTD)  solver  to  optimise  source  and  susceptor 
placement  inside  a  perfectly  conducting  structure,  build¬ 
ing  upon  initial  work  completed  in  [2j. 

Genetic  algorithms  are  briefly  introduced  in  section  2  and 
then  the  Finite-Difference  Time-Domain  method  is  intro¬ 
duced  in  section  3.  Section  4  describes  the  setup  of  the 
computer  simulation  and  section  5  discusses  the  results 
of  these  simulations,  finally  section  6  presents  conclusions 
from  the  work. 

2  Genetic  Algorithms 

GAs  are  a  stochastic  search  mechanism  with  thdr  opera¬ 
tion  firmly  rooted  in  natural  selection  and  survival  of  the 
fittest,  [3]  and  [4].  GAs  have  been  shown  to  provide  ro¬ 
bust  search  and  optimisation  in  complex  spaces,  [3]  and 
[5].  They  use  simple  operations  on  a  population  of  individ¬ 
uals,  which  lead  to  an  emergent  evolution  of  an  individual 
or  individuals.  Each  individual  in  the  population  repre¬ 
sents  a  potential  solution  to  the  given  problem  scenario 
and  as  such  is  evaluated.  After  an  individual  has  been 
evaluated  a  figure  of  merit  (FoM)  is  attributed  to  the  in¬ 
dividual.  This  FoM  is  a  measure  of  how  well  the  individual 
solves  the  problem.  GAs  can  lead  to  the  optimal  solution 
to  a  problem,  more  often  however,  they  are  used  to  op¬ 
timise  a  problem  towards  an  improved  solution  facilitat¬ 
ing  a  trade-off  between  excessive  computational  time  and 
meaningful  results.  The  GA  process  can  be  summarised 
graphically  as  shown  in  figure  1. 

The  GA  methodology  used  in  this  work  is  the  micro 
genetic  algorithm  (/^GA).  The  fiGA  technique  has  been 
shown  to  reach  the  optimal  area  of  multi-modal  solution 
spaces  earlier  than  conventional  GA  methods,  with  min- 
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Figure  1:  GA  Process  Flowchart. 

imal  computing  time,  [6].  It  has  already  been  demon¬ 
strated  that  the  pGA  can  be  successfully  linked  with 
an  FDTD  solver,  [7],  for  optimisation  in  EM  problems. 
Initially  in  a  GA  a  fixed  size  population  is  created  and 
populated  with  randomly  generated  individuals  of  a  fixed 
length,  this  length  is  determined  by  the  encoding  of  the 
problem  parameters.  The  population  size  for  the  pGA  is 
generally  kept  small  e.g.  five,  this  is  different  from  the 
classic  GA  where  the  population  size  is  typically  much 
larger,  in  the  order  of  100s.  These  individuals  are  all 
po.ssible  solutions  to  the  same  given  problem.  They  are 
asse.ssed  and  each  one  given  a  score,  the  FoM,  which  is 
then  used  to  determine  the  likelihood  of  reproduction  into 
the  next  generation.  The  FoM  is  analogous  to  fitness  in 
the  natural  world  where  fitter  individuals  survive  longer 
and  hence  have  a  better  chance  of  continuing  their  ge¬ 
netic  lineage.  The  particularly  successful,  higher  scoring 
individuals  may  be  reproduced  more  than  once  into  the 
next  generation.  In  the  /xGA  used  here  tournament  se¬ 
lection  is  employed,  this  is  perhaps  the  most  effective  for 
many  application  types  as  it  has  been  shown  to  provide 
better  convergence  towards  a  solution  in  the  initial  stages 
of  optimisation,  [8].  Tournament  selection  works  by  ran¬ 
domly  choosing  N  members  of  the  population  and  in  this 
instance  as  in  many  others  AT  =  2.  These  individuals  are 
then  pitched  against  each  other  to  determine  which  has  the 
better  FoM,  the  winner  of  the  tournament  being  selected 
to  be  a  member  of  the  new  population;  this  is  repeated  un¬ 
til  the  new  population  is  filled.  In  conjunction  with  tour¬ 
nament  selection  elitist  reproduction  is  also  employed,  this 
guarantees  that  the  best  individual  from  a  population  is 
present  in  the  next  population,  hence  preserving  the  best 
current  solution  to  the  given  problem  and  maintaining  a 
good  genetic  stock.  After  reproduction  the  individuals  un¬ 
dergo  crossover.  Here  two  randomly  picked  individuals  are 
mated  together,  swapping  information  between  their  chro¬ 
mosomes,  In  a  classic  GA  mutation  is  also  applied,  how¬ 
ever  in  the  f.iGA  mutation  is  turned  off.  Mutation  serves 
to  alter  the  genetic  makeup  of  these  new  offspring  with  a 
small  fixed  probability.  Once  again  mimicking  nature,  mu¬ 
tation  can  lead  to  either  a  detrimental  or  beneficial  effect 
on  performance.  Mutation  allows,  in  a  classic  GA,  random 
search  to  take  place  and  hopefully  leads  to  the  avoidance 


of  local  optima.  Once  the  population  has  converged  to  a 
determined  optimal  value  then  further  GA  operating  is  no 
longer  needed  and  the  GA  can  be  halted.  Due  to  the  small 
population  size  used  in  the  fiGA  premature  convergence  is 
often  seen  to  be  happening,  to  prevent  this  a  restart  mech¬ 
anism  is  used.  This  restart  mechanism  involves  checking 
for  convergence  in  the  current  generation  of  individuals 
and  then  restarting  the  next  generation  with  only  the  elite 
individual  and  new  randomly  created  individuals.  The  use 
of  a  restart  operator  ensures  random  search  takes  place 
and  leads  fiGA  away  from  local  optima.  Convergence  is 
checked  for  in  this  application  by  measuring  the  changes 
in  the  chromosomes,  the  individuals,  between  generations. 
Once  there  -are  few  changes  between  generations  then  the 
population  can  be  considered  to  be  converged,  the  setting 
of  this  limit  is  a  parameter  that  the  user  has  control  over. 
An  important  component  of  any  GA  code  is  the  device  by 
which  random  numbers  are  generated,  or  to  be  more  cor¬ 
rect  pseudo-random  numbers.  This  is  usually  a  portable 
pseudo-random  number  generator  (PRNG)  that  produces 
the  same  sequence  of  random  numbers  for  a  given  seed 
value.  The  quality  of  the  PRNG  is  an  important  factor. 
The  area  of  PRNG  research  is  vast  and  not  going  to  be,  for 
brevity,  discussed  here.  Common  references  on  the  subject 
are  [9],  [10],  [11]  and  [12]. 

A  clear  difference  can  be  seen  here  between  optimisation 
based  on  calculus  methods  and  GA  based  optimisation. 
GAs  make  use  of  interim  performance  in  the  optimisation 
problem,  calculus  based  methods  are  only  concerned  with 
optimisation  toward  an  overall  optimal  point  and  do  not 
“remember”  any  interim  performances. 

3  The  FDTD  Method 

The  FDTD  method,  [13],  [14],  has  become  one  of  the  most 
popular  methods  for  solving  Maxwell’s  equations.  It  is  a 
volumetric  domain  decomposition  technique  that  is  second 
order  accurate  in  space  and  time,  easily  implemented  in 
software  and  accurately  models  the  physical  world.  It  is  a 
widely  used  method  for  EMC/EMI  work,  [15],  as  well  as 
radar,  bioengineering  and  antenna  analysis.  The  FDTD 
method  is  described  widely  in  the  literature,  [16],  and  so 
only  a  very  brief  description  is  given  here. 

Maxwell’s  equations,  in  differentia]  form,  equations  (1)  - 
(2),  are  simply  modified  to  central-difference  equations, 
discretised,  and  implemented  in  software.  The  equations 
are  solved  in  a  leap-frog  manner,  [14],  Le.  tHe  electric  field 
is  solved  at  a  given  instant  in  tixpe,  then  the  magnetic  field 
is  solved  at  the  next  instant  in  time,  and  the  process  is 
repeated  over  and  over  again.  In  equations  (1)  and  (2)  H 
and  E  have  their  usual  meanings. 


11 

I 

(-VxB+^h), 

(1) 

il 

iv  X  H  -  -E. 
e  e 

(2) 

Using  a  three  dimensional  Cartesian  co-ordinate  system 
we  can  now  write  out  the  vector  components  of  the  curl 
operator.  Given  below  is  only  one  of  these  components, 
namely  the  x  component  of  the  H  field, 
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Figure  2:  Enclosure  geometry  shovnng  the  representations 
of  the  internal  structure  and  dielectric  slab,  all  sizes  are 
in  millimetres,  mm. 
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where  p  is  magnetic  permeability  and  p'  is  a  magnetic  loss 
parameter.  It  is  the  complete  set  of  coupled  partial  dif¬ 
ferential  equations,  six  in  total,  that  are  the  fundamental 
basis  of  the  numerical  FDTD  algorithm.  The  discretised 
form  of  equation  (3)  is  shown  below: 


+D,{m)  (  )  .  (4) 


where  n  is  the  time  step  under  consideration  and  i,3,k 
are  the  three  orthogonal  spatial  co-ordinates.  The  integer 
array  MEDIA  defines  material  conditions  for  each  field 
vector  component.  This  allows  the  medium  at  each  point 
to  be  mapped  out.  The  arrays  DayDb  are  magnetic  field 
update  coefficients,  there  are  corresponding  electric  field 
update  coefficients  also.  A  full  explanation  of  all  of  these 
equations  is  presented  in  [13]. 


4  Simulation  Setup 


4.1  Physical  Problem  Description 

The  problem  geometry  is  that  of  a  simplified  metal  box, 
which  has  perfectly  electrically  conducting  (PEC)  walls, 
an  internal  structure  and  DS  representation.  The  internal 
structure  is  modelled  as  two  PEC  sheets  in  the  top  corner 
of  the  geometry.  The  DS  can  be  thought  of  as  a  represen¬ 
tation  of  the  substrate  of  a  printed  circuit  board  (PCB). 
The  problem  geometry  is  shown  in  figure  2.  The  DS  has  a 
relative  permittivity,  Cr,  of  4  and  possesses  unity  magnetic 
permeability  pr. 


4.2  fxGA  and  FDTD  Setup 

The  /iGA  used  is  based  on  the  implementation  by  Carroll, 
[17],  modified  to  accommodate  the  given  task  of  moving 
the  source  point  and  susceptor  point  on  the  DS.  The  two 
points  are  specified  in  a  three  dimensional  Cartesian  sys¬ 
tem  with  the  y  ordinate  being  held  constant  as  this  rep¬ 
resents  the  surface  of  the  DS.  Binary  encoding  is  used  in 
the  pGA  and  this  leads  to  a  chromosome  string  length 
of  24  bits,  4  bits  per  a:,  y  and  z  co-ordinate  of  the  source 
and  susceptor.  It  should  be  noted  that  it  would  be  pos¬ 
sible  to  omit  the  y  ordinate  from  the  chromosome  how¬ 
ever  the  memory  saving  would  be  insignificant  especially 
when  compared  to  the  coding  effort  required  to  convert 
the  GA  for  this  one  specific  case.  Care  must  be  taken  to 
avoid  the  y  ordinate  being  altered  during  reproduction  and 
crossover,  this  is  achieved  using  a  uniform  crossover  oper¬ 
ator.  Uniform  crossover,  [18],  has  also  been  found  to  give 
faster  convergence  than  single  point  crossover  in  a  /iGA, 
[6],  [19].  The  population  size  is  mmntained  at  five.  These 
co-ordinate  values  are  passed  to  the  FDTD  solver  which 
computes  the  resulting  field  distribution  inside  the  enclo¬ 
sure.  The  peak  electric  field  at  the  susceptor  point  is  re¬ 
turned  as  the  FoM  to  the  pGA,  no  fitness  scaling  is  used  as 
tournament  selection  is  the  method  of  selection  employed 
in  this  GA  implementation,  [18].  The  pGA  attempts  to 
minimise  this  FoM,  i.e.  minimise  the  peak  electric  field  at 
the  susceptor  point. 


The  FDTD  solver  is  based  on  the  equations  given  in  [13]. 
The  DS  is  one  cell  thick  and  has  one  cell  of  free  space 
between  itself  and  the  enclosure  wall.  The  PEC  sheets 
are  modelled  as  infini^ly  thin.  A  soft  Er  dir^ted  sinu¬ 
soidal,  continuous  wave,  ideal  Hertzian  dipole  of  2.Q5GHz 
frequency  and  amplitude  of  10ym“^  is  u^ed  as  the  source 
in  the  FDTD  simulation.  The  domain  is  discretised  into 
dimensions  of  7.2727mm  in  the  x  and  z  directions  and 
7.2mm  in  the  y  direction.  Based  on  these  dimensions  and 
the  material  within  the  computational  domain  the  time 
step  size  for  the  solver  is  set  at  the  Courant  limit,  in  this 
instance  approximately  14ps. 


The  initial  results  are  obtained  from  using  the  G A  to  place 
a  source  and  susceptor  point  relative  to  each  other  on  the 
surface  of  a  dielectric  slab  (DS)  to  achieve  the  lowest  elec¬ 
tromagnetic  coupling  between  the  two.  The  complete  sim¬ 
ulation  is  implemented  in  Fortran  77  compiled  on  a  SUN 
ENTERPRISE  server  with  512MB  of  RAM  utilising  4  of 
8  processors. 


Initially  two  optimisation  cases  are  run;  a  non-lossy  and 
then  a  lossy  dielectric  slab,  in  the  lossy  case  the  conduc¬ 
tivity  is  set  to  0.04iSm”^.  After  this  a  multi-objective  sim¬ 
ulation  is  run  where  two  objectives  are  to  be  optimised, 
these  being  the  reduction  in  coupling  between  source  and 
susceptor  points  at  two  specific  frequencies. 
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5  Results 

5,1  Single  Objective  Optimisation 

The  progress  of  the  optimisation  process  can  be  seen  in 
figure  3  where  the  results  of  both  the  non-lossy  and  lossy 
simulations  are  presented.  This  figure  shows  the  evolution 
of  the  best  individual  in  the  simulation,  the  elite  indi¬ 
vidual.  While  this  individual  remains  the  best  the  FoM 
remains  at  the  same  value,  when  a  superior  one  is  bred,  a 
step  change  occurs  indicating  the  presence  of  a  new  elite 
value. 


Figure  3:  Generation  optimisation  of  both  the  non-lossy 
and  lossy  cases. 

The  final  value  for  the  non-lossy  case  is  0.0541  Vm“^  from 
an  initial  value  of  0.1172Vm“^  and  the  final  value  for  the 
lossy  case  is  0.0240Vm“^  from  a  value  of  0.0685Vm~^. 
Both  of  the  curves  are  normalised  to  the  value  of 
0.1172Vm“^  for  ease  of  comparison.  As  expected  the  val¬ 
ues  of  field  strength  for  the  lossy  case  are  lower  than  those 
for  the  non-lossy  case  as  in  the  lossy  case  energy  is  dissi¬ 
pated  in  the  dielectric  slab.  Figure  3  also  shows  us  that 
the  lossy  case  reaches  its  best  FoM  at  generation  66,  ear¬ 
lier  than  the  non-lossy  case  which  reaches  its  best  FoM 
at  generation  91.  This  is  a  random  dynamic  of  the  GA 
and  attributable  to  the  solution  surface  topology  and  the 
choice  of  PRNG  used  with  the  fiGA.  The  number  of  gen¬ 
erations  was  limited  to  100  as  experience,  with  this  code 
has  shown  that  there  is  a  minimal  return  on  computing 
time  by  exceeding  this  point.  Testimony  to  this  is  given 
by  the  fact  that  running  the  code  to  200  generations  gives 
only  a  final  value  of  0.0507Vm”^  for  the  non-lossy  case 
and  0.0204 Vm“'  for  the  lossy  case,  a  marginal  increase 
in  performance  for  a  doubling  of  computation  time.  This 
embraces  one  philosophy  of  using  a  GA,  namely  to  pro¬ 
duce  an  acceptable  improwment  toward  an  optimal  point 
for  a  limited  amount  of  effort. 

The  final  field  distribution  on  the  surface  of  the  DS  can 
be  seen  in  figure  4  for  the  non-lossy  case  and  in  figure  5 
for  the  lossy  case.  These  plots  are  generated  by  applying 
a  peak  hold  to  the  Ez  component  at  each  location  on  the 
surface  of  the  DS  throughout  the  simulation. 

The  final  positions  of  the  source  and  susceptor  points  on 


the  DS  are  indicated;  for  the  non-lossy  case  these  are  at 
(2,16)  and  (3,35)  respectively,  and  for  the  lossy  case  the  co¬ 
ordinates  are  (2,35)  and  (23,11).  These  co-ordinate  values 
given  are  of  course  on  the  x  and  z  axes  as  the  y  ordinate 
is  constant.  These  false  colour  plots  clearly  show  that  the 
fiGA  has  indeed  found  good  solutions  but  not  the  best 
solutions,  the  global  minimum  point  in  figure  4  is  at  (2,2) 
with  a  value  that  is  2AdB  down  on  the  position  found  by 
the  fiGA.  In  figure  5  the  global  minimum  is  at  (12,14)  with 
a  value  that  is  l.4dB  down  the  pGA  position.  The  final 
values  found  by  the  pGA  are  57.6dB  and  64.6dB  down  on 
the  global  maximum  for  the  non-lossy  and  lossy  cases  re¬ 
spectively.  It  should  be  noted  that  the  patterns  presented 
in  these  figures  are  ones  of  over  1  million  possible  patterns 
due  to  source  and  susceptor  placement,  finding  the  global 
minimum  on  each  of  them  would  require  direct  evaluation 
of  each  pattern,  the  fxGA  vastly  cuts  down  on  the  number 
of  evaluations  required  to  arrive  at  an  acceptable  result. 
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Figure  4:  False  colour  plot  of  final  field  distribution  on 
surface  of  dielectric  slab  after  non-lossy  simti/ofion.  The 
X  and  z  axes  are  marked  in  FDTD  cells.  Source  and  sus¬ 
ceptor  positions  are  indicated  on  the  figure. 

5.2  Multi- objective  Simulations 

Multi-objective  optimisation  (MO)  is  often  sought  in  prac¬ 
tice  as  often  compliance  in  one  particidar  objective  upsets 
the  compliance  of  another  objective.  MO  aims  to  achieve 
a  solution  that  can  not  be  improved  upon  .-simultaneously, 
in  each  of  the  objectives.  This  is  called  a  Pareto  optimal 
solution  and  the  set  of  all  these  solutions  is  called  the 
Pareto  optimal  set,  [20].  For  the  MO  simulations  some 
changes  had  to  be  made  to  the  GA-FDTD  code.  As  two 
objectives  were  being  optimised  a  method  of  measuring 
their  fitness  at  one  and  the  same  time  is  required,  to  ac¬ 
complish  this  the  method  of  objective  weighting  is  used, 
[20].  This  method  is  explained  mathematically  by  equa¬ 
tion  (5), 

Z  =  ^tt;i/i(x)  where  x  e  X.  (5) 

i=i 
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Figure  5:  False  colour  plot  of  final  field  distribution  on 
surface  of  dielectric  slab  after  lossy  simulation.  The  x  and 
z  axes  are  marked  in  FDTD  cells.  Source  and  susceptor 
positions  are  indicated  on  the  figure. 

In  equation  (5)  Z  is  a  scalar  valued  variable  that  is  a 
weighted  sum  of  the  individual  objectives,  where  there  are 
N  objectives  in  total.  The  function  /  is  the  function  that 
returns  the  FoM,  in  this  case  the  FDTD  solver,  and  x  rep¬ 
resents  the  parameters  of  the  function,  co-ordinate  values 
in  this  case.  The  feasible  region  of  solutions  is  denoted  by 
X.  The  sum  of  the  individual  weights  adds  to  1  and  each 
weight  lies  in  the  range  0  to  1  i.e.  0  <  toj  <  1.  It  is  the 
scalar  value  Z  that  is  optimised  by  the  GA.  The  method 
of  objective  weighting  is  an  easily  implementable  scheme 
that  produces  a  solution  from  the  Pareto  optimal  set. 

The  problem  setup  for  the  MO  optimisation  involves  the 
same  geometry  as  previously  with  the  only  change  being 
the  source.  A  Gaussian  pulse  modulated  onto  a  sine  wave 
is  now  used  as  the  source.  This  gives  a  symmetrical  spec¬ 
trum  about  the  carrier  frequency,  no  DC  component,  and 
a  bandwidth  controlled  by  the  length  of  the  pulse.  Its 
mathematical  form  given  in  [13]  is  described  by  equation 
(6), 

Eg  =  sin(27r/o(Ti  -  no) At).  (6) 

The  Gaussian  pulse  is  centred  around  frequency  /o  at  step 
no  with  a  1/e  characteristic  decay  of  ndecay  steps,  and  Af 
is  the  step  size.  The  time  series  and  frequency  response 
of  the  source  is  shown  in  figure  6.  The  centre  frequency 
for  the  source  is  chosen  at  1.500GHz  and  the  objectives 
for  optimisation  were  chosen  as  the  minimisation  of  the 
1.002GHz  and  2.004GHz  components  in  the  frequency 
response  which  is  computed  on  the  fly  as  the  FDTD  sim¬ 
ulation  progresses  as  detailed  in  [13].  These  frequencies 
were  chosen  as  they  tie  in  with  equipment  used  in  a  related 
measurement  study.  Only  non-lossy  simulations  were  run 
in  this  setup. 

The  resulting  frequency  response  from  the  Eg  component 
at  the  susceptor  point  the  simulations  can  be  seen  in  fig¬ 
ure  7.  The  dashed  line  in  figure  7  shows  the  initial  DFT 


before  the  optimisation  process  begins,  the  solid  line  in 
the  figure  displays  the  final  DFT  after  completion  of  the 
optimisation  process.  A  clear  difference  in  the  two  curves 
at  the  respective  objective  frequencies  can  be  seen.  Also, 
it  can  be  clearly  seen  that  the  field  values  at  the  objective 
frequencies,  on  the  final  DFT,  are  considerably  lower  than 
at  other  frequencies,  within  the  same  response,  excepting 
of  course  at  the  tails  of  the  response  where  there  is  mini¬ 
mal  energy  from  the  source.  The  values  of  the  objectives 
on  the  final  response  are  5idB  down  for  the  1.002GHz  ob¬ 
jective  on  the  max  value  in  the  response  and  the  value  at 
2.004GHz  is  44dB  down.  Relative  to  the  initial  DFT  the 
final  response  in  down  27dB  at  the  l,002GHz  objective 
and  down  ZSdB  at  the  2.004GHz  objective. 


Figure  6:  Time  series  of  the  pulsed  source^  upper  graph, 
and  its  DFT,  lower  gtaph.  This  is  the  source  characteristic 
used  in  the  multi-objective  optimisation  simulations. 


6  Conclusions 

A  method  of  using  a  pGA  in  conjunction  with  a  FDTD 
solver  to  facilitate  electromagnetic  optimisation  has  been 
shown.  The  conjoining  of  the  two  codes  allows  a  difficult 
design  problem  to  be  tackled,  namely  that  of  radiating 
and  susceptible  component  placement  within  an  arbitrary 
metallic- structure  to  minimise  coupling.  The  pGA  tech¬ 
nique  is  used  as  this  gives  a  significant  reduction  in  com¬ 
putation  time  with  little  loss  in  performance  of  the  final 
optimised  result.  Ensuring  that  the  ^GA  finds  the  global 
optimum  in  difficult  optimisation  problems  is  an  area  that 
needs  further  attention  but  the  ability  of  the  technique 
to  reach  a  near  optimal  solution  has  been  demonstrated. 
Both  single  ^nd  multiple  objectives  for  optimisation  have 
been  presented  with  promising  results  from  each.  It  is 
worth  noting  that  if  an  exhaustive  search  were  to  be  com¬ 
pleted  on  the  problems  presented  then  over  one  million 
simulations  would  be  required  taking  over  one  year  to  com¬ 
plete;  the  GA  based  search,  however,  took  approximately 
four  hours. 

The  design  optimisation  cases  presented  are  rather  sim¬ 
plistic  but  do  prove  the  concept  of  the  hybrid  code.  More 
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Figure  7:  DFTs  of  Ez  at  tiie  susceptor  point  of  the 
multi- objective  simulation  with  objective  of  minimising  the 
\.002GHz  and  2.004GHz  frequency  components  of  the  re¬ 
sponse.  Shown  are  the  initial  DPT  before  optimisation 
begins  (dashed  line)  and  the  final  DFT  after  optimisation 
(solid  line). 

challenging  problem  geometries  can  easily  be  accommo¬ 
dated  in  the  FDTD  code  with  little  burden  to  computing 
time  as  once  the  domain  has  been  discretised  it  is  the 
number  of  resulting  cells  that  chiefly  determines  the  com¬ 
putation  time.  It  is  envisaged  that  this  tool  can  be  used  to 
provide  design  rules  for  component  placement  within  an 
enclosed  resonant  cavity  and  that  it  can  can  also  be  used 
to  assess  specific  designs. 
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Abstract.  The  effect  of  discontinuities  at  edges 
and  at  feed-points  of  antennas  on  numerical 
convergence  rates  Is  investigated.  In  the  case 
of  edges,  higher  order  representations  of  the 
edge-mode,  expressed  with  the  aid  of  Hermite 
splines,  are  shown  to  provide  improved 
convergence  in  both  global  and  local  measures. 
When  using  a  magnetic  frill  to  excite  an  antenna, 
it  is  shown  that  when  the  current  representation 
allows  for  the  "charge  jump"  across  the  frill,  then 
convergence  Is  accelerated.  The  use  of  both 
sub-domain  and  entire-domain  functions  is 
explored. 

Introduction. 

The  treatment  of  edge  conditions  in 
electromagnetic  problems  is  generally  thought  to 
be  well  understood,  particularly  after  the 
publication  by  Meixner  [1]  of  the  relationship 
between  'edge  angle’  and  the  power  of  the 
distance  from  the  discontinuity.  An  extensive 
review  by  Van  Bladel  [2]  brought  together  much 
of  the  evidence  in  support  of  the  phenomenon’s 
existence.  Nevertheless,  the  model  proposed 
by  Meixner.  while  widely  accepted,  has  not 
benefited  from  further  development.  More 
recently,  Shen,  et  al.,  [3][4]  reported  on  an 
analytical  solution  for  the  current  and  charge 
distributions  at  the  ends  of  a  tubular  dipole. 
Their  work  showed  that  the  Meixner  method 
applied  only  at  the  very  ends  of  such  a  dipole. 
While  extending  our  knowledge  in  this  area,  th^r 
work  is,  nevertheless,  restricted  to 
current/charge  distributions  on  a  tubular  dipole. 

A  different  form  of  discontinuity  arises  at  the 
feed-point  of  an  antenna.  This  fact  is  less  well 
recognized.  The  author  knows  of  only  one- 
report  that  discusses  it  quantitatively.  In  this 
report,  Popovic,  Dragovic  and  Djordjevic  [5,  p38] 
pointed  out  that  “as  a  consequence  of  the 
annular  magnetic-current  frill  excitation,  the 
current  derivative  has  a  discontinuity  at  the  frill 
location”.  They  then  went  on  to  quantify  the 
charge  jump  that  occurs  in  this  situation. 


Unfortunately,  it  is  not  clear,  from  the  balance  of 
their  writing,  how  they  consistently  incorporated 
this  observation  into  their  work  in  a  general 
manner. 

Recently,  the  author  published  reports  [6][7]  that 
examined  the  use  of  higher  order  basis  functions 
in  modeling  a  tubular  dipole  antenna/scatterer. 
The  results  presented  in  these  reports  did  not 
support  the  expectation  that  higher-order  basis 
functions  would  provide  faster  convergence 
properties.  The  models  used  in  these 
investigations  employed  splines  up  to  fifth  order. 
They  also  incorporated  the  Meixner  edge 
condition.  Subsequent  analysis  has  shown  that 
the  application  of  the  simple  form  of  the  edge 
condition  actually  inhibited  the  convergence 
properties  of  the  higher  order  models.  The 
simple  form  of  the  Meixner  edge  condition  has 
only  one  degree  of  freedom  -  it's  amplitude.  In 
order  for  the  splines  covering  the  rest  of  the 
antenna/scatterer  to  meet  their  potential,  the 
number  of  degrees  of  freedom  in  the  edge 
model  were  increased  (in  this  study)  to  match 
those  of  the  adjacent  spline. 

As  part  of  the  analysis  referred  to  in  the 
preceding  paragraph  an  appreciation  for  the 
problem  associated  with  the  charge  jump  across 
the  magnetic  frill  developed.  The  problem 
exhibits  itself  through  a  slow  convergence 
and/or  oscillation  in  the_  input  admittance  of  the 
..antenna  as  the  model  size  is  increased.  ”  The 
slow  convergence  problem  was  found  to  occur 
for  both  entire-domain  and  sub-domain  basis 
functions. 

The  Basis  Functions. 

The  sub-domain  basis  functions  used  on  the 
main  portion  of  the  antenna/scatterer  are  all 
splines  and  are  shown  in  Table  I.  They  are  the 
same  as  those  used  in  [6][7].  They  include  the 
linear  spline,  which  acts  as  a  reference  for  the 
performance  of  the  higher-order  functions.  As 
the  results  when  using  the  3-term-sin/cos 
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Expansion  Function 

Mathematical  Description 

3-terni  sin/cos 

/(y)  =  4  +  5,sin(A:(>»  -  jy,))  +  q(cos(A:Cy  -  y,)  -  1) 
where  ^  jy  < 

Linear  Spline 

liy)  =  a,  +  b^iy  -  y,)  y,  <  y  < 

Quadratic  Spline 

liy)  =  a,  +  b,{y  -  y^)  +  q(y  -  y,f  y,  ^  y  ^  jy,., 

Cubic  Spline 

J^(y)  =  «;  +  bi(y  -  y,)  +  c^(y  -  y^f  +  d,{y  -  y^f 
where  y,  ^  y  ^  y^ 

Quartic  Spline 

Ky)  ^  a,  +  b.,{y  -  y,)  +  c^iy  -  y^f  +  4(p  -  y,f  +  e^y  -  yi)* 
where  y,  ^  y  <  y^ 

Table  1.  Definitions  ol 

F  the  splines  considered  in  this  study. 

N  N-\  N  N  ^ 

On*\  =  +  i,^A,  +  a^  (la)  =  2YcA  +  (1b) 

1=  fsl  1*1  /=! 

functions  were  almost  identical  to  those  found  function  can  be  performed  no  more  than  once 
when  using  the  quadratic  spline,  only  the  latter  while  still  maintaining  a  function  that  can  be 


are  reported.  The  reader  should  be  aware 
however  that  the  3-term-sin/cos  functions  were 
studied.  When  using  splines  such  as  those 
found  in  Table  I,  one  asserts  as  much  continuity 
at  the  junction  of  contiguous  cells  as  the 
functions  allow.  For  example,  a  quadratic 
spline,  which  is  of  order  3  and  degree  p  -  2, 

provides  C*  continuity.  Hence,  in  the  case  of  a 
quadratic  spline  both  the  amplitude  and  the  first 
derivative  can  be  matched.  By  this  means  the 
number  of  variables  is  reduced.  In  the  case  of 
the  quadratic  spline,  for  N  cells,  one  needs  to 
determine  only  the  free  variables 

N 

amplitude  and  first 

/=! 

derivative  at  one  extremity  are  related  to  those 
at  the  other  extremity  through  the  relationships 
given  in  equations  (1a)  and  (1b)  above. 

The  sub-domain  edge/end  basis  functions  used 
in  this  study  are  based  on  the  form  given  in  [2, 

p120]: 

E,  =  d"{a^  +  Oyd  +  + . )  (2) 

where  is  the  electric  field  normal  to  the 
edge/end,  d  is  the  distance  from  the  edge/end, 
V  is  a  number  greater  than  -1 .0  and  the  series 

Oo.flj.Oj.Oj, . are  constants  to  be  determined. 

Expression  (2)  needs  to  be  matched  up,  at  its 
cell  boundary,  with  each  of  the  functions  shown 
in  Table  I.  Thus  derivatives  of  this  function  are 
required.  Straightforward  differentiation  of  the 


integrated,  and  thus  might  be  employed  with  the 
quadratic  spline.  However,  it  would  fail  for 
higher  order  splines.  An  alternative,  and 
numerically  superior,  method  is  to  apply  the 
concepts  found  in  the  development  of  Hermite 
splines. 

Hermite  splines  are  different  from  most  splines 
found  in  computational  electromagnetic  studies 
which  are  concerned  only  with  amplitudes.  The 
splines  in  Table  I  are  examples  of  this  latter 
situation.  The  constants  at  the  open  ends  of 
these  splines  are  the  source  of  much  writing  on . 
ways  to  handle  them.  -  see,  for  example,  [8]. 
The  derivatives  of  these  splines  are  functions  of 
foese  constants.  In  the  development  of  Hermite 
splines  the  issue  of  derivatives  is  tackled 
directly.  Hermite  splines  are  typically  of  odd' 
degree,  including  p  =\,  and  are  defined  in 

terms  of  their  amplitude  and  (p  - 1)/2 
derivatives  at  each  end.  The  process  for 
developing  Hermite  cubic  splines  is  covered  in 
[9]  and  the  more  general  case,  with  applications 
to  CEM,  is  described  in  [10].  In  the  case  of 
equation  (2)  only  one  end  is  open  and  the 
development  for  a  two  term  expression  will  be 
described  next.  Two  terms  will  provide 

amplitude,  f(1),  and  derivative, 

information  that  can  be  used  directly  with  the 
quadratic  spline  described  earlier. 
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Equation  (2)  is  restated  as 
/>(«)  =  /(])//,(«/)  +  - 

(A  -  A)  5  z'  ^  A ,  A  is  the  width  of  the  cell  and 
h  is  identified  with  a  celi  boundary.  //(<  )  is  a 
polynomiai  of  order  2  selected  so  that  H^(\)  =  1 
and  =  0.  Similarly,  H^iu)  is  a  simiiar 

polynomial  selected  so  that  AfjO)  =  0  and 
Hi  ’(1)  =  1 .  When  these  conditions  are  applied, 
the  solution  is  given  by: 

P{u)  =  u'\(\  +  v)  -  yu)f(l) 

+w^(l  -  u)AfX\) 

This  provides  us  with  what  we  are  looking  for  -  a 
function  that  separates  the  amplitude  and  the 
first  derivative  at  the  open  end  and  one  that  can 
be  integrated  in  a  straightforward  manner.  This 
concept  can  be  extended  to  functions  of  higher 
orders,  and  is  done  in  this  study  for  use  with  the 
cubic  and  quartic  splines  of  Table  I. 

When  investigating  the  discontinuity  at  a 
magnetic  frill,  a  structure  that  possesses  no 
other  discontinuities  is  useful  for  investigative 
purposes.  A  circular  loop  is  such  a  structure 
and  is  examined  here.  The  ‘natural’  basis 
functions  to  use  in  this  situation  are  Fourier 
series.  Due  to  symmetry,  orily  cosine  terms  are 
typically  used.  The  expansion  series  is  then: 

N 

I  =  cos(/^)  0  <  <  2/r  (4a) 

1*0 

This  expansion  has  been  In  use  for  decades  and 
the  results  of  its  use  are  well  documented  [11]. 
When  one  examines  the  reported  results,  the 
lack  of  finality  In  the  convergence  of  the 
admittance  at  the  feed-point  is  notable.  Zhou 
and  Smith  [12]  attributed  this  poor  convergence 
to  the  use  of  a  delta-function  source.  Their 
results  using  a  magnetic  frill  source 
demonstrated  that  indeed  the  frill  source 
performed  better  than  the  delta-function  source. 
However  the  admittance  curves  still  did  not  fully 
converge  until  N  in  (4a)  was  extremely  large. 
This  is  surprising  given  the  ’natural’  suitability  of 


the  functions  in  (4a).  The  work  reported  here 
shows  that  the  model  in  (4a)  is  deficient  for  use 
when  dealing  with  the  magnetic  frill  as  it  does 
not  account  for  the  ‘charge  jump*  across  the  frill. 
The  derivatives  of  (4a)  are  all  continuous,  and 
what  is  needed  is  a  set  of  functions  that  allows 
for  this  ‘charge  jump’.  Instead  of  discarding 
(4a),  one  can  modify  it  as  in  (4b)  below. 

I  =  a,  +  +  A,sin((2/  -  1)^/2))  (4b) 

i=\ 

The  derivatives  of  this  series  are  such  that  a 
'charge  jump*  can  be  modeled  properly.  As 
discussed  in  [10],  in  the  immediate  vicinity  of  the 
frill  it  is  the  odd  derivatives  of  the  current  that 
must  be  In  anti-phase  on  opposite  sides  of  the 
frill. 

Solution  Tools 

The  equations  used  in  the  analysis  that  follows 
are  based  on  the  Electric  Field  Integral 
Equation,  EFIE,  formulation.  In  particular,  in  the 
case  of  the  linear  dipole,  the  equation  is  Hallen’s 
derivation  [13],  shown  In  (5)  below.  In  (5)  /  is 
the  desired  current,  G  is  the  Green’s  function, 

R  =  ^(2-zy -f(2asin(Y))^ ,  the  radius  is  a ,  C 

is  a  constant  to  be  determined  and  Ei  is  the 
incident  excitation.  The  dipole,  which  is  an 
open  circular  cylinder  with  an  infinitesimally  thin 
wall  thickness,  Is  located  with  its  center  at  the 
origin  of  a  cylindrical  coordinate  system  with  its 
length  aligned  with  the  z  axis. 

When  investigating  the  loop,  the  derivation  due 
to  Mei  [14],  specialized  for  a  loop  is  used  -  see 
(6)  below.  The  center  of  the  loop  is  located  at 
the  origin  of  g  cylindrical  coordinate  system,  with 
s  -  Rf,^.  The  loop  itself  lies  in  the  plane 
with  its  feed-point  at  ^  =  0 .  hi  this  case,  the 
expression  for  R  in  G  is 
R  =  +  2RI(\  -  cos^^)  +  2aRf,(l  -  cos^')cos^‘ 

,  where  Is  a  local  variable  around  the  circular 
cross-section  of  the  loop  element.  Here  the 
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Extended  Boundary  Condition  is  used  and  the 
net  electric  field  on  the  axis  of  the  loop  is  zero. 
The  numerical  solutions  for  I  in  equations  (5) 
arid  (6),  which  are  expressed  in  the  usual 
ZI  =  V  system,  are  obtained  with  the  aid  of  the 
Boundary  Residual  Method,  BRM,  [15][16]. 
This  is  a  Least  Squares  Method  variant  that  has 
two  important  elements.  The  first  significant 
element  is  that  the  entries  in  the  Z  and  V 
matrices  are  calculated  as  they  would  be  for  a 
point-matching  method.  However  the  location 
of  the  matching  points  is  not  arbitrary. 
Specifically  they  are  located  at  the  nodes  of  an 
appropriate  integration  rule  such  as  Gauss- 
Legendre.  Furthermore,  the  number  of  sample 
points  exceeds  the  number  of  free  variables, 
typically  by  a  ratio  of  2:1  or  higher.  This  means 
the  Z  matrix  is  rectangular  and  the  matrix 
equations  are  solved  using  QR  or  SVD  methods 
[17].  The  second  important  element  is  that  the 
rows  in  the  Z  and  V  matrices  are  weighted  by 
the  weights  from  the  same  integration  rule  used 
to  establish  the  nodes  mentioned  earlier.  This 
has  the  important  result  of  minimizing  the 
residuals  integrated  over  the  entire  surface 
being  studied.  Formally,  this  statement  is  min 

pI^  =  [Z/  -  .  This  is  an  important 

measure.  It  reports  on  the  performance  of  the 
solution  over  the  complete  surface  and  hence  it 
is  referred  to  here  as  a  global  measure.  For 
comparative  purposes,  it  is  best  to  normalize 
this  function.  The  normalized  residual  is: 


The  Dipole  as  a  Scatterer. 

In  order  to  calculate  the  currents  bn  a  straight 
dipole  excited  by  a  plane  wave  equation  (5)  is 
evaluated  in  conjunction  with  the  basis  functions 
shown  in  Table  I.  The  boundary  condition, 
I(±h)  =  0 ,  is  satisfied  in  one  of  three  ways,  in 
the  first  method  a^  and  are  set  to  zero.  In 
the  second  approach,  the  two  end  cells  support 
the  basic/conventional  end  function 

/(z')  =  a^^(h  -  and  then  enforce 

=  a^  =  af,^^ .  The  derivatives,  b,  and 
and  higher,  are  left  unconstrained.  In  the  third 
approach,  the  special  Hermite  end  splines 
defined  earlier  are  used  and  all  applicable 
constraints  are  enforced.  In  the  case  of  the 
linear  spline,  the  second  and  third  approaches 


utilize  the  same  function.  For  the  quadratic 
case,  a  two  term  end  spline  is  used,  for  the 
cubic  case  a  three  term  end  spline  is  used  and 
the  quartic  spline  employs  a  four  term  end 
spline.  The  plane  wave  is  polarized  parallel  to 
the  axis  of  the  dipole  and  is  normally  incident. 

Figure  la  shows  the  resuits  for  the  current  at  the 
center  of  a  half-wave  dipole,  with  a  radius  of 
0.007A .  As  the  order  of  the  spline  is  increased 
the  results  without  the  inclusion  of  any  special 
end  sections  show  only  a  small  benefit  from  the 
inciusion  of  higher  order  terms.  The  inclusion  of 
the  conventional  end  spline  provides  an 
improvement,  particularly  for  the  lower  order 
splines,  interestingly,  the  convergence  cun/es 
for  all  four  splines  look  very  similar.  It  is  only 
when  the  Hermite  type  end  splines  are  added 
that  the  real  benefit  of  using  high  order  splines  is 
observed. 

Figure  1b  reports  the  results  for  the  Normalized 
Residual.  It  is  observed  that  the  conventional 
end  function  improves  the  performance  of  the 
lower  order  basis  functions.  However,  the 
performance  of  the  Hermite  end  functions  is 
consistently  better. 

The  Loop  Antenna. 

As  mentioned  above,  the  circular  loop  is  used  to 
investigate  the  discontinuity  associated  with  a 
magnetic  frill.  The  basis  functions  defined  in 
equations  (4a)  and  (4b)  are  used  in  conjunction 
with  equation  (6)  to  evaluate  the  currents  at  the 
feed-point  of  the  antenna.  When  using 
equation  (4b)  the  value  of  N  is  only  half  that  of 
the  corresponding  value  used  with  equation 
(4a).  The  circumference  of  the  loop  is  one 
wavelength  and  n  =  21n(2/rf?4/<7)  =  10.0, 
where  is  the  radius  of  the  loop  and  a  is  the 
Tadius  otthe  loop  element.  .  '  ' 

The  convergence  curves  for  the  central  current 
are  shown  in  Figure  2a  and  the  associated 
normalized  residuals  are  shown  in  Figure  2b. 
The  results  show  clearly  the  superiority  of  a 
series  that  accommodates  the  'charge  jump’  at 
the  feed-point,  thereby  confirming  the  presence 
of  the  phenomenon.  Only  32  terms  are 
considered  in  these  graphs.  If  the  number  of 
terms  were  increased  substantially,  the  results 
for  the  cosine  series  only  would  eventually 
converge  to  those  observed  with  the  use  of  the 
supplemented  basis  functions. 


Log10(Normalized  Residual)  |l|  .  milliamps 
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Figure  1b.  Plots  of  the  normalized  residual  on  a  half-wave  dipole  excited  by  a  plane  wave. 
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7  =  ^0+  ^aiCos(mz'lh)^h  ~  \z'\  (^3) 

i=I 

I  -  Oq  +  ^(flf.cos(/;r2y2A)  +  bfSm({2i  -  I);r|z'|/2/?)yA  -  \z]  (8b) 

/=i 


The  Dipole  as  an  Antenna.  Discussion. 


The  dipole  antenna  incorporates  both  of  the 
types  of  discontinuity  already  discussed.  Here, 
it  has  radius,  a,  of  0.007/1 .  The  frill  has 
dimension  of  f)/a=2.3  where  b  is  the  outer  radius 
of  the  frill.  The  combination  of  these  two 
discontinuities  is  explored  initially  using  sub- 
domain  functions.  The  sub-domain  functions 
are  the  same  as  those  used  earlier  with  a  dipole 
scatterer.  However,  at  the  center  of  the  dipole, 
the  phase  In  the  terms  representing  the  first  and 
third  derivatives  is  shifted  180”  when  these 
derivatives  exist.  The  impact  on  the 
convergence  of  the  current  at  the  feed-point  can 
be  seen  In  the  results  reported  in  Figure  3a. 
The  results  for  the  associated  normalized 
residuals  are  shown  In  Figure  3b.  It  is  clear, 
from  both  sets  of  data,  that  the  phase  reversal  of 
the  odd  order  derivatives  at  the  frill  contributes 
significantly  to  an  improved  model  for  the 
current. 

Entire-domain  functions  have  long  been  used  to 
model  the  current  on  a  linear  dipole  antenna.  A 
study  that  compared  ten  such  functions  was 
published  recently  [18].  That  study  concluded 
that  entire-domain  functions  performed  poorly  on 
the  linear  dipole.  Historically,  however,  the 
choice  for  an  entire-domain  function  on  a  linear 
dipole  has  been  one  of  various  forms  of  a 
Fourier  series.  In  the  absence  of  a  better 
alternative,  such  a  series  Is  used  here.  The 
particular  form  of  the  Fourier  series  used  is 
shown  In  (8a)  above.  In  order  to  accommodate 
the  “charge  jump"  at  the  magnetic  frill,  a 
modified  form  of  the  latter  was  investigated  - 
this  is  shown  in  (8b)  above. 

The  results,  for  the  input  admittance  of  a  half¬ 
wavelength  antenna,  using  these  two  series  are 
shown  in  Figure  4a.  The  corresponding 
findings  for  the  normalized  residuals  are  shown 
in  Figure  4b.  When  using  equation  (8b)  the 
value  of  N  is  only  half  that  of  the  corresponding 
value  used  with  equation  (8a),  Both  figures 
clearly  demonstrate  the  utility  of  incorporating 
terms  that  address  the  Issue  of  discontinuity 
across  the  frill. 


End  Effects.  When  considering  the  results 
reported  in  Figure  la,  one  must  bear  in  mind 
that,  for  the  linear  spline,  the  conventional  end 
term  and  the  Hermite  end  term  are  one  and  the 
same.  Adding  a  single  end  term  to  any  of  the 
splines  produces  an  observable  improvement  in 
the  convergence  of  the  current  at  the  center  of 
the  dipole  compared  with  not  using  one  at  all. 
This  is  particularly  true  for  the  lower  order 
splines.  As  mentioned  earlier,  the  convergence 
curves  obtained  when  using  a  single 
conventional  end  term  all  appear  to  have  the 
same  convergence  characteristics.  This 
suggests  that  the  end  term  is  limiting  the 
convergence  process.  When  the  higher  order 
end  terms  are  included,  to  the  maximum  order 
possible,  the  convergence  process  is  speeded 
up  significantly.  When  viewed  from  a  global 
perspective,  as  reported  in  Figure  1b,  the 
conventional  end  term  is  important  in  the  case  of 
the  linear  spline.  For  the  higher  order  splines 
the  benefits  become  less  as  the  order  of  the 
spline  is  increased.  However,  when  the  higher 
order  end  terms  are  included,  the  global  results 
are  distinctly  better  for  all  splines. 

Magnetic  Frill.  The  importance  of  properly 
Incorporating  “charge  jump”  properties  Into 
models  of  current  flowing  on  an  antenna  excited 
by  a  magnetic  frill  is  dramatically  illustrated  in 
Figures  2-4.  The  significant  oscillations  seen 
when  employing  a  Fourier  series  containing  only 
cosine  terms  is  an  indication  of  Its  lack  ot 
suitability  -  for  both  the  Idbp  and  the  linear 
dipole!  To  be  sure  that  the  oscillations  were  the 
result  of  the  feed  mechanism  and  nothing  else, 
the  excitation  was  changed  from  the  magnetic 
frill  to  a  plane  wave.  The  oscillations 
disappeared  in  the  case  of  the  loop  and  were 
highly  attenuated  on  the  linear  dipole.  This 
allowed  the  conclusion  that  the  oscillations  are 
indeed  a  result  of  using  the  magnetic  frill 
excitation  and  not  something  else.  The  addition 
of  the  sine  terms  in  the  current  model,  while 
Improving  the  convergence  results  significantly, 
does  raise  the  condition  number  of  the 
impedance  matrix  from  the  region  of  ~  10^  up  to 
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Figure  3a.  Plots  of  admittance  on  a  dipole  excited  by  a  magnetic  frill  for  three  splines. 
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Figure  4a.  Plots  of  the  admittance  of  a  half-wave  dipole  excited  by  a  magnetic  frill 


Figure  4b.  Plots  of  the  normalized  residual  on  a  half-wave  dipole  excited  by  a  magnetic  frill 
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~  10'°  when  more  than  12  terms  are  used  in  the 
series. 

The  iinear  dipoie  problem  was  also  solved  using 
a  Method-of-Moments  solution,  using  the 
standard  Pocklington  equation.  This  was  done 
as  a  check  on  the  testing  method  and  the 
integral  equation  formulation.  The  testing 
functions  were  pulses  of  equal  dimensions,  the 
basis  functions  were  again  the  Fourier  series  of 
(8a)  and  (8b)  and  the  Pocklington  equation  was 
regularized  per  the  procedure  described  in  [19]. 
The  results  of  using  (8b)  converged  significantly 
faster  than  when  using  (8a)  while  the  oscillations 
noted  in  connection  with  Rgure  4  were  almost 
completely  attenuated.  These  findings  provided 
further  confirmation  of  the  need  to  accommodate 
the  discontinuity  at  the  feed-point. 

The  presence  of  a  high  condition  number  in  the 
entire  domain  solutions  that  include  a  series 
expansion  of  the  absolute  value  of  sine  terms 
makes  them  unattractive.  Consequently,  the 
effect  of  incorporating  a  single  term  only  was 
investigated.  In  the  case  of  the  loop  this  term 
was  sin(|(^'|/2) .  In  the  case  of  the  dipole  it  was 

simply  Iz*! .  The  associated  results  are  shown 

in  Figures  2  and  4  respectively.  Use  of  just 
these  single  terms  contributes  significantly  to 
alleviating  the  convergence  problem.  The 
condition  numbers  were  comparable  to  those 
observed  when  using  the  cosine  terms  only. 
This  suggests  that  one  judiciously  selected  term 
could  significantly  speed  up  the  convergence 
rate  in  these  situations.  However,  the  ones 
tested  here  are  not  those  functions  -  otherwise 
they  would  exhibit  normalized  residual  curves  at 
least  as  steep  as  those  shown  by  the  addition  of 
the  sine  series. 

The  basic  EFIE  for  a  loop  is  shown  in  equation  - 
(9a).  After  some  manipulation,  this  assumes 
the  form  shown  in  (9b).  When  the  current  is 
represented  by  (4a)  the  form  in  (9c)  results. 
The  last  term  on  the  right-hand  side  of  this  latter 


equation  evaluates  to  zero.  But  this  is  the  term 
that  represents  the  “charge  jump".  Therefore 
the  representation  of  (4a)  is  inadequate.  A 
term,  or  terms,  that  permit(s)  the  current 
derivative  at  the  origin  to  be  discontinuous  is 
necessary.  The  inadequacy  of  the 
representations  of  (4a)  and  (8a),  is  believed  to 
be  the  main  source  of  the  oscillations  seen  in 
the  relevant  convergence  curves. 

In  1980,  Richmond  [20j  reported  that  the 
incorporation  of  a  term  that  represented  an 
outgoing  wave  on  an  infinite  dipole  improved 
convergence  considerably  on  a  finite  dipole. 
This  was  confirmed  in  [18J,  where  it  was  also 
noted  that  the  condition  numbers  were  high. 
This  outgoing  term  is  essentially  a  log  term  at 
the  feed-point  and  consequently  presents  a 
discontinuity  there.  It  is  possible  that  this 
discontinuity  had  more  to  do  with  the 
improvement  realized  by  Richmond  than  the 
form  of  the  term  itself. 

Conclusions. 

This  study  was  motivated  by  the  apparent  failure 
of  high-order  basis  functions  to  improve 
convergence  rates  when  compared  with  basis 
functions  of  lower  order  [6][7j.  As  a  result  of 
this  work,  it  is  believed  that  the  reasons  for  this 
failure  originate  with  inadequate  current 
modeling  in  the  vicinity  of  discontinuities. 
Specifically,  it  has  been  demonstrated  that  to 
achieve  superior  convergence  properties  one 
must  bear  in  mind  three  things. 

1.  Current  models  in  the  vicinity  of 
discontinuities  must  support  the  same  level  of 
continuity  as  that  supported  by  the  models 
employed  away  from  the  discontinuities. 

2.  When  an  isolated  excitation  source, 
such  as  a  magnetic  frill,  is  part  of  the  system, 
the  current  model(s)  must  have ‘sufficient 
flexibility  to  accommodate  issues  such  as 
"charge  jump”  in  the  excitation  region. 

3.  The  above  comments  apply  to  both  sub- 
domain  and  entire-domain  models. 
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Abstract:  This  paper  describes  a  three-dimensional 
electromagnetic  propagation  model  for  signal 
power  prediction  in  a  non  line  of  sight  urban  area 
located  in  Singapore.  The  model,  which  is 
implemented  using  the  Ohio  State  NEC-BSC  V42 
basic  scattering  code  and  is  based  on  ray  theory 
and  uniform  theory  of  diffraction  (UTD)  takes  into 
account  first-order,  second  order,  and  third  order 
effects  including  triple  reflections  and  diffractions. 
The  simplified  propagation  model  of  Raffles  Place, 
the  central  business  district  of  Singapore,  uses 
more  than  360  geometrical  structures  and  is 
compared  with  data  measured  at  937.6  MHz  along 
a  drive  route  at  the  site.  The  results  from  the 
propagation  model  produced  reasonable  agreement 
with  the  measurements  showing  that  a  simplified 
model  using  UTD  can  be  used  to  simulate  the  gross 
features  of  electromagnetic  scatter  in  an  urban  area. 
It  is  shown  that  use  of  penetrable  dielectric 
building  walls  are  also  necessary  for  accurate 
prediction  of  radio  wave  propagation  in  urban  areas 
under  non  line  of  sight  conditions.  It  is  also  found 
that  third  order  scattering  effects  can  be  dominant 
in  non  line  of  sight  situations  and  may  be  necessary 
for  accurate  predictions. 


1  INTRODUCnON 

As  cellular  mobile  communications  systems 
become  more  widely  used,  there  is  an  increasing 
need  to  develop  reliable  planning  tools  for  base 
station  deployment.  Part  of  this  growth  is  due  to 
frequency  reuse  whereby  the  same  frequencies_are 
utilized  in  geographically  distinct  areas  or  c3ts 
such  that  signal  path  loss  over  distance  can  be 
>  exploited.  Frequency  reuse  provides  for 
significantly  greater  capacity.  One  of  the 
drawbacks  with  frequency  reuse  however,  is  the 
potential  for  interference  between  nearby  cells.  As 
cells  become  smaller  in  size,  the  use  of  reliable 
planning  tools  becomes  even  more  necessary.  In 
these  environments,  accurate  signal  strength 
prediction  plays  an  important  role  in  controlling 
intercell  interference  and  providing  coverage  for 
low-power  handheld  units.  The  use  of  reliable 
planning  tools  means  that  the  use  of  expensive  and 
time-consuming  measurements  in  the  field  can  be 


minimized.  These  tools  can  also  provide  a  cost 
effective  means  of  anticipating  problems  at  the 
early  stages  of  base  station  deployment  and  of 
optimizing  parameters  such  as  antenna  location  and 
orientation. 

Several  techniques  have  been  developed  [1-6] 
for  analysis  of  multipath  and  other  propagation 
effects  in  macrocellular  and  microcellular 
environments.  In  several  of  these  approaches,  ray 
tracing  or  launching  techniques  are  combined  with 
GTDAJTD  (Geometrical  Theory  of  Diffraction 
AJniform  Theory  of  Diffraction)  to  analytically 
determine  the  sum  of  all  incident,  reflected,  and 
diffracted  fields  at  the  receiver  location.  The  use 
of  GTD/UTD  is  appropriate,  since  the  frequencies 
used  in  mobile  communications  (e.g.  >  900  MHz) 
are  usually  high  enough  such  that  the  assumption 
that  the  features  of  most  building  structures  are 
much  larger  than  the  wavelength  used  becomes 
valid.  In  [2],  a  computer  based  propagation  model 
is  used  to  obtain  a  time  delay  representation  of  the 
received  signal.  In  [3,  4],  a  multislit  waveguide 
with  randomly  distriWed  gaps  between  the  sides 
of  buildings  is  considered  as  a  model  of  straight 
streets  under  line  of  sight  conditions.  In  [5],  a  ray 
launching  technique  is  used  to  detect  the 
intersection  of  a  ray  with  an  object  and  the 
diffracted  and  reflected  rays  are  then  computed. 
This  process  is  repeated  until  a  maximum  number 
of  reflections  is  exceeded.  In  [6],  2D  FDTD  is 
corrected  to  account  for  spherical  wave  spreading 
and  applied  to  a  simple  outdoor  environment.  Each 
building  in  the  model  consists  of  a  perfectly 
_  conducting  core  with  a  dielectric  coating. 

In  this  paper,  a  particular  urban  area  was 
selected,  a  three-dimensional  model  was  created 
and  the  propagation  characteristics  were  simulated 
using  NEC-BSC  V4.2  (Numerical  Electromagnetic 
Code-Basic  Scattering  Code)  [7],  The  model 
presented  includes  contributions  to  the  received 
signal  from  all  possible  propagation  paths, 
including  ground  and  wall  reflections  from 
diffracted  and  specularly  reflected  signals  both  in 
the  LOS  (Line  of  Sight)  and  NLOS  (Non  Line  of 
Sight)  regions  using  GTD/UTD.  This  code  is 
structured  such  that  all  of  one  type  of  scattered 
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field  or  interaction  is  computed  at  one  time  to 
increase  efficiency.  Buildings  having  an  arbitraiy 
layout  can  be  modeled  even  if  the  streets  are  not 
rectilinear  and  the  buildings  walls  on  each  side  of 
the  streets  are  not  coplanar. 

The  purpose  of  this  research  is  to  investigate 
how  accurately  UTD  could  predict  electromagnetic 
scattering  in  a  non  line  of  sight  urban  environment 
using  a  simplified  model  consisting  of  plates  and 
curved  surfaces.  In  order  to  evaluate  the 
limitations  of  the  model,  predictions  from  the 
model  considering  f,  and  3'^  order  scattering 
mechanisms  were  compared  to  CW  test  data 
measurements  from  a  private  cellular  provider. 
This  paper  is  organized  as  follows:  Section  2 
describes  the  urban  propagation  model,  Section  3 
compares  the  measured  vs.  predicted  results,  and 
Section  4  summarizes  the  results  and  provides 
some  directions  for  future  research. 

2  URBAN  PROPAGATION  MODEL 

The  goal  of  this  investigation  was  to  determine 
how  well  UTD  could  be  used  to  predict 
electromagnetic  scattering  in  a  heavily  urbanized 
environment  using  a  simplified  model  of  the  site. 
Raffles  Place,  the  central  business  district  of 
Singapore  was  chosen  as  an  appropriate  site  for 
this  study  due  to  its  high  building  density.  Other 
reasons  for  choosing  Raffles  Place  were: 

1.  The  building  layout  does  not  lie  completely  in 
a  rectangular  grid.  This  allows  for 
investigations  to  be  carried  out  in  an  area 
without  a  rectangular  grid  layout  and  having  a 
somewhat  random  building  orientation. 

2.  Drive  routes  can  be  investigated  which  arc 
mainly  non-line-of-sight  (NLOS). 

3.  The  building  heights  in  this  area  range  from  5 
to  282  meters.  The  effects  of  diffraction  from 
the  edges  of  rooftops  and  sides  of  buildings 
can  be  investigated. 

Figure  1  shows  a  three-dimensional  view  of 
the  geographical  layout  of  Raffles  Place  used  in  the 
UTD  model  with  the  drive  route  shown.  The  drive 
route  begins  along  Cecil  Street  and  follows  onto 
Collyer  Quay  up  to  Fullerton  Square  then  turns  left 
onto  Battery  Road  where  it  ends  on  Chulia  Street. 
The  area  covered  is  about  500  m  by  500  m  and  the 
length  of  the  drive  route  is  1.12  km.  The  data 
predicted  along  this  route  is  compared  with  the 
measured  drive  data  collected  by  MobileOne,  a 
private  cellular  network  provider  in  Singapore.  In 
all  field  measurements,  GPS  was  used  to  record  the 


location  of  the  measured  signal  power  along  the 
route.  The  following  input  model  parameters  are 
required  in  the  NEC-BSC  UTD  code  [7]: 

1 .  Coordinates  of  the  plates,  curved  surfaces,  and 
other  structures  necessary  to  simulate  building 
exterior  and  roofs. 

2.  Estimate  of  the  electrical  permittivity, 
conductivity,  and  thickness  of  the  building 
exteriors. 

3 .  Transmitter  and  receiver  locations. 

4.  Transmitter  and  receiver  antenna  patterns. 

5.  Operating  frequency  used. 

6.  Types  and  order  of  UTD  interactions  to  be 
included  in  the  simulation. 

The  coordinates  of  the  building  exteriors  were 
extracted  from  an  aerial  map  by  using  an  overlay 
having  a  fixed  origin  selected  on  the  map.  The 
comers  of  the  buildings  were  measured  relative  to 
this  origin.  Most  building  shapes  were  modeled 
using  only  plates,  which  were  assumed  to  be  flat 
surfaces  having  average  electrical  permittivity  and 
conductivity.  Structures  having  low  curvature 
sections  were  modeled  using  multiple  plates  to 
form  a  piecewise  linear  approximation  of  that 
section.  When  the  radius  of  curvature  of  the 
structure  was  substantial,  cylinders  or  truncated 
cones  were  used.  Elliptic  cylinders  were  also  used 
to  model  some  building  structures. 

The  ground  plane  was  simulated  using  a  single 
homogeneous  half  space  of  a  relative  dielectric 
permittivity  Er  =  10  and  conductivity  a  =  0.001 
S/m,  which  can  be  considered  reasonable  for  a 
typical  ground.  Preliminary  studies  suggest  that 
the  predicted  signal  power  near  the  ground  is 
somewhat  insensitive  to  the  ground  dielectric 
constant  and  conductivity.  The  seawater  surface 
was  also  included  in  the  model  and  was  modeled  as 
a  dielectric  coated  plate  with  a  relative  permittivity 
of  Er=  80  and-electrical  conductivity  o  =  4  S/m 

The  model  database  used  in  the  UTD 
propagation  model  consists  of  360  plates  and  5 
cylinders  to  represent  the  64  buildings  and  other 
structures  located  at  Raffles  place.  The  number  of 
plates  in  the  model  could  easily  exceed  500  if  the 
entire  area  is  to  be  accurately  modeled.  However, 
following  the  assumption  that  effects  from 
structures  like  trees,  lamp  posts,  telephone  boxes. 
Mass  Rapid  Transit  station  buildings,  etc  can  be 
ignored,  these  structure  were  not  included  in  the 
model  [8].  Every  attempt  was  made  to  define  the 
model  in  more  detail  in  areas  where  it  was 
illuminated  strongest  by  the  antenna.  Because  of 
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the  wide  variation  in  building  heights,  the  detailed 
description  in  the  model  is  not  limited  to  the  area 
local  to  the  transmitter.  Since  the  model  used 
omnidirectional  antennas  for  the  transmitter  and 
receiver,  the  entire  area  including  but  not  limited  to 
the  region  between  the  transmit  antenna  and  the 
receiver  field  point  was  modeled.  As  a  result  of  the 
above,  the  number  of  plates  required  for  modeling 
was  quite  large  to  account  for  all  the  possible  ray 
paths. 

The  building  exteriors  were  modeled  as 
dielectric  plates  to  approximate  materials,  such  as 
glass  windows,  absorbing  panels,  concrete,  etc. 
Research  is  still  underway  to  produce  useful 
approximations  for  the  thickness,  permittivities  and 
conductivities  of  these  delectric  layers  [8],  The 
choice  of  a  good  approximation  for  the  electrical 
properties  is  of  utmost  importance,  since  the 
Fresnel  reflection  coefficients  depend  on  these 
values.  These  coefficients  dictate  the  reflection 
coefficient  of  the  materials  influencing  the 
modeled  results. 

The  buildings  located  at  Raffles  Place  were 
constructed  using  a  variety  of  materials,  such  as 
concrete,  glass,  and  steel.  It  is  difficult  to  identify 
the  electromagnetic  properties  of  each  individual 
building  exterior  and  it  is  assumed  that  all  the 
building  walls  have  electrical  property  values 
averaged  over  those  of  glass  and  concrete.  The 
building  walls  were  modeled  using  plates  and 
cylinders  coated  with  dielectric  on  one  side  to 
simulate  the  building  exterior  and  a  perfect  electric 
conductor  on  the  other.  The  material  properties 
were  chosen  as  Ei-  =  15,  a  =  7  S/m,  and  thickness  = 
0.5  meter  [9,  10].  The  accuracy  of  the  model  is 
limited  by  the  lack  of  accurate  building  and  street 
databases  and  the  approximation  of  building 
exteriors  to  be  “smooth”  flat  surfaces  having  an 
average  relative  permittivity  and  conductivity. 

~  Both  the  transmitter  and  receiver  antennas 
were  simulated  using  vertically  polarized  half¬ 
wave  dipoles.  The  transmitter  input  power  was 
32.4  dBm.  The  transmitter  dipole  antenna  was 
placed  3  m  above  the  machine  room  on  the  rooftop 
of  Ocean  Building,  a  building  of  about  120  m  in 
height.  The  receiver  was  placed  2  m  above  the 
ground  to  simulate  mounting  of  a  vertical 
monopole  on  top  of  a  car.  The  measured  data  was 
taken  at  intervals  of  between  28  -  59  m  using  a 
GPS  based  measurement  system.  The  simulated 
data  is  sampled  in  20  m  intervals  with  and  without 
spatial  averaging.  All  UTD  scattering  mechanisms 
were  selected  during  the  simulation.  The  NEC- 


BSC  UTD  code  however,  allows  for  the  individual 
scattering  terms  to  be  selected  to  investigate  their 
behavior  along  the  route. 

3  COMPARISON  OF  THE  MEASURED 

VS.  PREDICTED  RESULTS 

Three  sets  of  simulation  cases  were 
considered.  The  first  simulation  case  considers 
only  the  1®‘  and  2"^*  order  interaction  terms  listed  in 
Table  1 .  The  second  simulation  case  considers  up 
to  and  including  the  order  interactions.  The 
additional  3**  order  interaction  terms  are  listed  in 
Table  2.  For  these  first  two  cases,  the  building 
walls  only  reflect  and  diffract  the  radio  waves 
using  the  material  dielectric  properties  given.  The 
third  simulation  case  considers  f‘,  7^*,  and  3^ 
order  interactions  listed  in  Tables  1  and  2  using 
penetrable  dielectric  plates  as  building  walls.  For 
this  case,  radio  waves  are  both  reflected  and 
transmitted  through  the  plates  to  simulate 
propagation  through  building  structures. 

Figure  2  shows  a  comparison  between  the 
predicted  signal  powers  resulting  from  the  first  two 
simulation  cases  (1®‘  and  2"**  order  interactions 
only,  and  1®\  and  3^**  order  interactions)  and  the 
measured  fields.  The  f  *  and  2?*^  order  interaction 
terms  listed  in  Table  1  result  in  a  total  of  5015  ray 
paths  along  the  drive  route.  Table  2  lists  the  3^ 
order  interaction  terms  resulting  in  an  additional 
42850  ray  paths.  The  curve  showing  all  1st, 
and  3**  order  interactions  appears  to  give  better 
agreement  with  the  measured  field  data. 

The  greatest  discrepancy  between  the 
measured  and  predicted  data  for  cases  1  and  2 
occurs  over  the  0.2  km  stretch  located  immediately 
before  point  D  near  Fullerton  Square  (Figure  1). 
The  predicted  results  for  the  case  of  up  to  3"*  order 
interactions  underestimate  the  measured  data  by 
more  than  20  dB.  This  is  difficult  to  explain  using 
the  model  because  there  are  no  nearby  structures 
located  in  this  area  on  which  the  rays  can  get 
scattered.  Along  the  drive  route  to  the  right  of 
Collyer  Quay  is  an  open  area  immediately  adjacent 
to  Marina  Bay.  One  possible  explanation  for  the 
discrepancy  between  the  measured  and  predicted 
data  is  the  lack  of  radio  wave  propagation  through 
buildings  since  modem  office-building  interiors  are 
composed  of  mostly  air.  Studies  have  shown  that 
radio  waves  may  penetrate  into  buildings  with  only 
between  9-11  dB  attenuation  at  900  MHz  [11]. 
Our  model  was  modified  to  allow  the  building 
walls  to  be  penetrable  by  radio  waves  such  that  the 
total  attenuation  through  the  building  structures  is 
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between  12-15  dB.  The  properties  of  the  building 
walls  were  changed  to  =  15,  0.0023  S/m,  and 

thickness  -  1  m  and  were  simulated  as  penetrable 
dielectric  plates.  Shown  in  Figure  2  is  the 
measured  vs.  predicted  signal  power  along  the 
drive  route  in  Figure  I  using  penetrable  walls.  The 
correspondence  between  measured  and  predicted 
power  is  much  better  than  for  the  cases  where  the 
building  walls  only  provide  reflection  and 
diffraction.  It  can  be  observed  that  use  of 
penetrable  building  walls  results  in  II  dB  less 
variation  in  the  predicted  field  along  the  drive  route 
than  for  the  case  of  3"^*  order  scatter  only.  Table  3 
lists  the  RMS  error  and  standard  deviation  for  each 
simulation  case. 

Additional  observations  can  be  made  from  the 
measured  vs.  predicted  signal  power  results  shown 
in  Figure  2.  At  location  A  on  Cecil  Street,  the  peak 
in  the  measured  received  field  may  be  attributed  to 
3*^^^  order  interactions  since  this  peak  cannot  be  seen 
in  the  2"*^  order  results.  For  instance,  in  the 
prediction,  a  3^^*  order  reflection  occurs  between  the 
Singapore  Land  Tower  building,  the  OUB  Center, 
and  a  second  reflection  from  the  Singapore  Land 
Tower  (Figure  3).  Multiple  interactions  between 
the  Ocean  Building  and  the  Ocean  Towers  also 
occur  before  the  ray  arrives  at  the  receiver  location. 

At  location  B,  the  measured  signal  power  is  at 
the  lowest  level  along  the  drive  route  and  the 
number  of  2"**  and  order  ray  paths  appears  to  be 
minimum.  At  location  B  the  only  contribution  to 
the  signal  power  at  the  receiver  location  appears  to 
be  from  the  diffraction  from  the  rooftop  of  the 
Ocean  building.  This  location  is  also  the  point 
closest  to  the  transmitter  along  the  drive  route. 

At  location  C,  the  predicted  results  that  include 
up  to  2"^  order  interactions  only  indicate  a  drop  of 
nearly  50  dB  in  the  received  signal  strength  due  to 
obscuration  by  an  overpass  located  along  the  route. 
For  the  predicted  results  that  include  3^  order 
interactions,  the  effect  of  the  overpass  is  minimal. 
The  measured  signal  power  also  indicates 
essentially  no  effect  due  to  the  overpass .  The 
effects  of  the  3"*^  order  interactions  appear  to  be 
dominant  directly  under  and  immediately  beyond 
the  overpass.  Multiple  rays  that  originate  at  the 
transmitter  location  may  be  diffracted  from  the 
rooftop  edge  on  top  of  the  Ocean  Building  and  then 
diffract  and  reflect  from  the  north  side  of  the  John 
Hancock  building.  These  rays  are  then  reflected 
from  the  sides  of  both  the  Union  Overseas 
Shopping  Center  and  a  restaurant  and  then  finally 
arc  reflected  from  the  side  of  the  Hitachi  Center 


(see  Figure  4)  to  the  receiver  location  immediately 
beyond  the  overpass  location. 

Probably  the  most  interesting  observation  in 
the  measured  and  predicted  data  is  shown  at 
location  D.  Here,  the  measured  and  predicted 
signal  powers  are  at  their  maximum  value  over  the 
entire  route.  Location  D  is  in  front  of  the 
Singapore  Land  Tower  on  Battery  Road.  The 
received  signal  power  is  about  12  dB  above  the 
average  signal  power  along  the  entire  route.  At 
this  location,  the  transmit  antenna  on  the  Ocean 
Building  is  obstructed  by  several  buildings.  The 
received  signal  appears  to  come  from  several  ray 
paths  that  first  begin  as  reflections  from  the  OUB 
Center  and  UOB  Plaza  and  then  are  reflected  from 
the  Singapore  Land  Tower  and  Standard  Chartered 
buildings  to  the  receiver  location.  This  area  of 
Battery  Road  appears  to  behave  like  a  parallel-plate 
waveguide  so  that  the  rays  from  several  locations 
are  guided  into  a  relatively  narrow  street  bounded 
by  the  Singapore  Land  Tower  on  one  side  and 
Standard  Chartered  Bank  and  May  Bank  on  the 
other  side  (see  Figure  5). 

At  the  end  of  the  route  and  furthest  from  the 
transmitter  on  Battery  street  (location  E)  is  the 
location  of  the  largest  discrepancy  between  the  2"^ 
and  3^**  order  predicted  data.  The  measured  signal 
power  in  the  area  around  location  E  however,  is 
nearly  constant.  The  predicted  data  however, 
shows  almost  a  50  dB  variation  for  the  1**  and  2"*^ 
order  interactions  and  nearly  27  dB  of  variation  for 
the  case  of  up  to  3*^*^  order  interactions.  The 
predicted  data  shows  a  2"*^  order  interaction  that  is 
the  Esult  of  a  specular  reflection  from  the  City 
House  building  and  then  diffraction  between  the 
edges  of  the  Exchange  Building  and  Golden  Shoe 
Car  Park  (Figure  6)  to  the  receiver  location.  The 
measured  power  increase  is  only  about  1  dB  at  this 
location.  Rays  reflected  from  the  Republic  Plaza 
and  the  UOB  Plaza  channel  into  this  area  resulting 
in  a  peak  in  the  received  signal  power.  Like  the 
results  shown  at  location  D,  this  area  may  behave 
like  a  parallel-plate  waveguide  where  multiple 
reflections  occur  between  UOB  Plaza  on  one  side 
and  the  Sinsov  building  on  the  other  side  across 
Battery  Street. 

4  CONCLUSION 

In  this  paper,  a  simplified  model  using  UTD 
with  ray  tracing  and  consisting  of  360  plates  and  5 
cylinders  was  used  to  predict  high-frequency 
electromagnetic  wave  propagation  in  an  urban 
environment.  Taking  into  account  single,  double, 
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and  third  order  UTD  terms  interacting  between 
building  walls  and  comers,  the  gross  features  of  the 
electromagnetic  response  were  accurately 
predicted.-  The  three-dimensional  model  was 
developed  using  arbitrary  building  layouts  modeled 
using  simple  stmctures  constructed  using  plates  to 
represent  flat  surfaces  and  cylinders  for  curved 
surfaces. 

Several  sets  of  UTD  simulations  were 
performed  and  the  resulting  predicted  signal 
powers  were  compared  with  measurements  along 
the  drive  route  located  in  Raffles  Place,  the  central 
business  district  of  Singapore.  Comparisons 
between  the  predictions  and  measurements  suggest  * 
that: 

1.  When  there  is  no  line  of  sight  (NLOS), 
diffraction  cannot  be  neglected.  The  effects  of 
diffraction  from  the  edges  of  the  rooftop 
appeared  to  be  significant  at  the  transmitter 
location.  and  order  interactions  alone 
may  not  be  sufficient  for  accurate  predictions. 
Inclusion  of  higher  order  scattering  is 
important  for  the  accurate  prediction  of  the 
scattering  in  NLOS  situations.  Inclusion  of  3^^ 
order  interactions  also  result  in  9  dB  less 
variation  than  2"**  order  interactions  in  the 
predicted  signal  power  along  the  drive  route. 

2.  For  NLOS  cases,  the  fields  at  the  receiver 
location  often  result  from  scattering  on 
structures  local  to  the  receiver  location. 
These  fields  may  also  result  from  reflections 
from  large  structures  or  surfaces  far  away 
from  the  receiver  location.  This  can  result  in 
propagation  paths  that  may  not  be  intuitively 
obvious.  These  higher  order  scattering  terms 
(3'*^  order  and  above)  may  require  a  larger 
model  to  account  for  all  possible  ray  paths. 

3.  Lack  of  sufficient  structures  in  the 
projiagation  model  and  inaccuracies  in 
antenna  positions,  route  locations,  may 
introduce  errors  into  the  predicted  signal 
powers.  Another  source  of  potential  error  is 
scattering  from  cars,  trucks,  or  other  vehicles 
located  near  the  receiver  and  not  included  in 
the  model.  It  was  also  assumed  that  the 
dielectric  constant  and  conductivity  of  the 
building  exteriors  throughout  the  model  was 
constant.  These  simplifications  introduce 
extra  error. 

4.  UTD  can  also  produce  artifacts  in  the 
predicted  fields.  At  some  locations 


considering  up  to  order  interactions  only 
result  in  artifacts  not  observed  in  the  measured 
data.  At  those  locations,  the  3”*  order  scatter  is 
dominant  over  the  f  and  order  scatter 
effects. 

5.  At  two  locations  along  the  drive  route,  it  was 

found  that  the  peak  in  the  Eceived  signal 
power  might  be  the  result  of  a  parallel-plate 
waveguide  effect  caused  by  multiple 

reflections  between  adjacent  buildings 

separated  by  a  street  or  road. 

6.  Effects  of  radio  wave  propagation  through 
building  walls  and  structures  are  necessary  in 
resolving  differences  between  measured  and 
predicted  data  under  NLOS  conditions.  Use  of 
penetrable  building  walls  in  the  model  reduces 
the  RMS  error  in  the  prediction  by  13  dB. 

Future  work  may  include  incoiporation  of 
higher  order  (e.g.  4*’’  order)  interaction  terms  in  the 
propagation  study.  This  may  be  useful  for  the 
prediction  of  the  signal  power  near  Fullerton 
Square.  The  dielectric  properties  of  the  building 
exteriors  are  also  very  important  and  may  have  a 
large  effect  on  the  magnitudes  of  reflections  from 
buildings.  Future  work  includes  additional 
modeling  of  propagation  through  buildings. 
Analysis  of  these  effects  will  be  helpful  in  , 
resolving  differences  between  measured  and 
predicted  signal  power  data.  Also  of  interest  are 
the  cross-polarized  fields  scattered  in  non  line  of 
sight  locations,  for  the  study  of  polarization 
diversity  antenna  schemes. 
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Table  2. 3^*^  order  interactions  considered  in  case  2. 


RPRPRP 

Triple  reflection  from  plates 

RPRPDP 

Double  Plate  Reflection  -  Edge  Diffraction 

RPRPVP 

Double  Plate  Reflection  -  Comer  Diffraction 

RPDPRP 

Plate  Reflection  -  Edge  Diffraction  -  Plate 
Reflection 

RPVPRP 

Plate  Reflection  -  Comer  Diffraction  -  Plate 
Reflection 

DPRPRP 

Edge  Diffraction  -  Double  Plate  Reflection 

Measured  vs.  Predicted  Signal  Power 


Figure  2.  Measured  and  predicted  signal  power  along  drive  route. 


Table  3.  RMS  error  and  standard  deviation  between  measured  and  predicted  power  along  drive  route. 


Case 

Order  of  interactions 
included  in  model 

RMS  Error 

Standard  Deviation 

1 

js.  2"<I 

33.0  dB 

13.2dB 

2 

2nd  3rd 

22.5  dB 

13.1  dB 

3 

2"*^,  3^*^  with  penetrable 
walls 

9.2  dB 

8.1  dB 

Ellis:  Wireless  Propagation  in  Non  Line  of  Sight  Urban  Areas 
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Standard  Chartered 
Bank 


UOB  Plaza 


Figure  6. 


Major  Ray  contribution  at  Location  E. 
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Abstract  -  In  the  framework  of  photonic  ctystal’s  band 
structure  calculations,  we  present  a  novel  way  -  based  on 
several  advanced  techniques  for  searching  and  tracing 
eigenvalues  with  the  multiple  multipole  program  -  to  compute 
these  diagrams  automatically,  efficiently,  and  with  a  high 
accuracy.  Finally,  we  validate  the  results  for  some  well  known 
test  cases. 

Keywords  -  band  diagram;  multiple  multipole  program;  photonic 
crystal;  eigenvalue  analysis. 

I.  Introduction 

HOTONIC  Crystals  (PhCs)  were  proposed  in  1987  by  E. 
Yablonovitch  [1]  at  the  University  of  California,  as  an 
optical  counterpart  to  semiconductors,  i.e.,  PhCs  should 
provide  a  photonic  bandgap  in  the  same  way  that  a 
semiconductor  possesses  an  electronic  bandgap.  In  fact,  PhCs 
are  rarely  found  in  nature.  Exceptions  are  opals  and  butterfly 
wings.  However,  thanks  to  nano-technology  it  has  become 
possible  to  fabricate  artificial  PhCs  in  the  last  decade.  These 
PhCs  essentially  consist  of  a  periodic  assembly  of  dielectric 
scatterers,  i.e,,  there  is  a  strong  link  to  the  well-known 
structures  of  grating  theory.  One  of  the  important  differences 
between  PhCs  and  semiconductors  is  the  size  of  the  unit  cell. 
For  a  semiconductor,  one  has  a  3D  grid  consisting  of  identical 
atoms,  i.e.,  the  lattice  constant  in  all  three  directions  of  the 
crystal  is  in  the  order  of  the  diameter  of  an  atom,  whereas  the 
cell  size  of  a  PhC  is  in  the  wder  of  half  an  optical 
wavelength,  i.e.,  much  larger.  From  this  fact,  one  expects  that 
the  Photonic  Integrated  Circuits  (PICs)  based  on  PhC  concept 
must  be  much  larger  than  the  traditional  semiconductor  ICs. 
But  -  because  of  the  macroscopic  size  of  the  PhC*s  unit  cell  - 
one  has  much  more  freedom  in  introducing  and  fabricating 
defects  in  a  PhC  than  in  a  semiconductor.  Note  that  a 
semiconductor  becomes  interesting  from  the  technological 
point  of  view  only  when  a  few  impurities  or  defects  are 
introduced  and  the  same  holds  for  PhCs.  Doping  traditional 
semiconductors  is  a  rather  statistical  way  of  introducing 
defect  atoms  in  a  semiconductor  and  therefore,  the  building 
blocks  of  semiconductor  are  relatively  large  blocks  of  material 
with  a  specific  doping.  These  blocks  obviously  consist  of 
many  atoms.  When  designing  a  "doped"  PhC,  one  can 
precisely  position  and  fabricate  all  defects  with  a  high  degree 


of  freedom  -  at  least  it  is  expected  that  this  can  be  done  in  the 
near  future  [2]. 

Although  the  variety  of  PhC  structures  that  might  be 
fabricated  one  day  seems  to  be  almost  infinite  and  although 
many  interesting  structures  were  already  proposed  (various 
types  of  waveguides,  sharp  bends  in  waveguides  without  any 
reflection,  couplers,  resonators,  etc.)  or  even  fabricated  on  a 
prototype  level,  one  currently  cannot  say  what  kind  of  PhC 
structures  will  be  favored.  At  the  moment,  one  can  neither 
know  the  materials  that  are  best  suited  for  PhCs  -  it  is  well 
known  that  a  large  dielectric  contrast  is  required  for  obtaining 
a  bandgap,  which  somehow  limits  the  materials  that  may  be 
used,  but  there  is  no  unique  choice  at  all  -  nor  what  kind  of 
geometry  {2D  crystals  or  3D  crystals  [3-5],  symmetry,  shape 
of  the  scatterers)  is  most  appropriate.  Thus,  there  is  a  strong 
need  for  theoretical  investigations  and  simulations  of 
potential  structures.  The  first  step  of  such  investigations 
consists  in  the  computation  of  the  band  diagrams  of  perfect 
PhCs  without  any  defects.  The  goal  is  to  find  structures  that 
may  easily  be  fabricated  and  exhibit  a  broad  band  gap,  i.e.,  a 
frequency  range  where  no  electromagnetic  waves  are  allowed 
to  propagate  within  the  ciystal.  In  order  to  find  the  band  gap, 
one  must  compute  the  band  diagram  of  the  lowest  order 
modes  of  the  PhC.  This  is  essentially  an  eigenvalue  problem 
that  exhibits  several  special  cases  that  may  cause  difficult 
numerical  problems,  especially  when  one  is  designing  a 
procedure  for  the  aut^atic,  efficient,  accurate;  and  reliable 
computation  of  the*^complete  band  diagrams  for  arbitrary 
structures. 

Currently,  the  most  frequently  used  approach  is  the  Plane 
Wave  Method  (PWM)  that  mainly  approximates  the 
electromagnetic  field  by  a  superposition  of  plane  waves  [6- 
10].  It  is  well  known,  that  this  method  has  a  problematic 
convergence  [11-13,  10].  Other  methods  that  were  used  for 
PhCs  are  the  Finite  Difference  Time  Domain  (FDTD) 
[14,15],  the  transfer  matrix  method  [16],  the  Finite  Element 
Method  (FEM)  [17],  and  the  Boundary  Element  Method 
(BEM)  [18].  In  the  following,  we  apply  the  latest  version  of 
the  Multiple  Multipole  Program  (MMP)  [19]  implemiented  in 
the  MaX-1  software  [20]. 
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In  order  to  obtain  efficient,  reliable,  and  accurate  results,  we 
carefully  analyze  the  numerical  problems  that  may  occur  and 
introduce  several  new  techniques.  For  reasons  of  simplicity 
we  focus  on  the  2D  case. 

The  remainder  of  the  paper  is  organized  as  follows:  A  commonly 
used  representation  of  PhCs  in  terms  of  their  band  diagram  is 
elucidated  in  Section  11.  In  Section  III  we  briefly  explain  the 
core  of  our  photonic  crystal  calculations,  when  MMP  is 
considered.  The  proper  framework  of  the  eigenvalue  search  is 
reported  in  Section  IV,  whereas  in  Section  V  a  successful 
automation  of  such  search  procedure  is  proposed.  A 
validation  of  our  band  structure  calculation  by  means  of 
various  test  examples  is  given  in  Section  VI.  And  finally,  we 
conclude  this  contribution  with  a  short  summary  in  Section 
Vll. 


II.  Definition  of  the  Band  Diagram 

As  an  introductory  example  let  us  consider  the  simple  case 
of  a  2D  PhC  consisting  of  dielectric  rods  arranged  in  a  square 
lattice  and  embedded  in,  e.g.,  air.  For  periodic  structures  it  is 
possible  to  apply  some  fundamental  theorems  from  solid  state 
physics.  The  original  lattice  for  this  crystal  is  given  on  the  left 
hand  side  of  Fig.  1 .  For  the  dielectric  constant  we  can  write 
e{r)  =  €(f  +  R)  where  R  is  one  of  the  original  lattice  vectors. 
According  to  Bloch’s  theorem  [6,  7]  for  the  modal  field  inside 
the  crystal  we  write: 

E=Ei^(?)  =  Ui^(r)e*’  ’  (,) 

=  +  ’  (2) 

Note  that  (1)  holds  not  only  for  the  electric  but  also  for  the 
magnetic  field.  Bloch’s  tlieorcm  may  be  proven  in  classical 
electro-dynamics  [6].  Important  consequences  of  this  theorem  are 
[6,7]: 

1.  c'*^  =  l,  i.c.,  k'R  =  N27r^  where  N  is  an  integer  -  the 
wave  vector  space  (reciprocal  space)  is  discrete, 

2-  E-  -  (/“)»  Ic.,  the  reciprocal  space  is 

periodic.  G  is  one  of  the  reciprocal  lattice  vectors. 

This  allows  us  to  define  the  so-called  reciprocal  lattice  "Space, 
spanned  by  the  reciprocal  lattice  vectors.  We  first  define  the 
original  lattice  vectors  as  follows: 

R^71A+7j,e,+vA  (3) 

where  e, ,  Cj*  are  three  independent  lattice  vectors  and 
77j,  7/3  arc  integer  numbers.  Note  that  ej  is  missing  in  2D 
crystals.  Similarly,  we  write  for  the  primitive  reciprocal 
lattice  vectors: 


^  “  ^\f\  ^2/2 


(4) 


If  we  want  to  construct  the  reciprocal  lattice,  we  can  use  [7]: 


7,=2;r 


e^  •  (fij  X  Cj) 


/3  =  2^- 


These  equations  are  derived  from  the  definition  of  the 
reciprocal  lattice  vector  space.  For  2D  crystals  (cylindrical 
geometry),  the  vector is  omitted  and  the  vector  is  the 


unit  vector  e,  along  the  cylinder  axis. 

From  the  equations  above,  we  can  conclude  that  the  discrete 
translational  symmetry  of  a  photonic  crystal  leads  to  the  fact 
that  modes  with  the  wave  vector  k  and  modes  with  the  wave 


vector  k  +  G  are  identical,  i.e.,  we  have  periodicity  also  in  the 
reciprocal  space.  A  special  representation  of  the  primitive  cell 
for  this  periodicity  is  called  the  first  Brillouin  zone  (l*‘  BZ).  It 
can  be  defined  as  a  zone  around  any  lattice  point  in  the 
reciprocal  space  with  points  that  arc  closer  to  this  lattice  point 
than  to  any  other  lattice  point. 


0 

? 

0 

0 

( 

y- 

0 

0 

0 

Fig.  1 :  The  original  (left)  and  reciprocal  (right)  lattice  for  a  2D  photonic  crystal 
(square  lattice).  Construction  details  for  reciprocal  lattice  are  given  in  the  text. 


Fig.  2:  C:onstruction  of  the  1”  Brillouin  zone  (solid  square),  its  irreducible  part 
(triangle  P-X-M)  and  characteristic  points  for  band  structure  computation  (P, 
X,  and  M). 

The  Brillouin  zone  construction  (using  Bragg’s  planes  - 
dashed  lines)  for  the  square  lattice  is  shown  in  Fig.  2. 
Because  of  the  high  degree  of  symmetry,  wc  need  to  analyze 
only  a  small  part  of  the  1’*  BZ.  This  part  is  called  the 
irreducible  BZ  (IBZ),  [6,  7],  In  the  case  of  periodic  structures. 
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Fig.  3:  The  unit  cell  of  the  photonic  crystal  with  dielectric  rods  arranged  on  a 
square  lattice. 

it  is  sufficient  to  perform  the  modal  field  analysis  in  the  area 
of  the  EBZ.  As  illustrated  in  Fig.  2  the  IBZ  for  a  square  lattice 
is  a  triangle  with  the  corners  F,  X,  and  M.  Since  the  maxima 
and  minima  of  the  eigenvalues  {resonance  frequencies)  are 
supposed  to  be  on  the  borders  of  the  IBZ,  it  is  sufficient  to 
trace  the  eigenvalues  along  the  sides  of  the  IBZ  in  order  to  get 
the  photonic  bandgaps.  Therefore,  the  standard  band  diagram 
consists  of  three  sections:  F-X,  X-M,  and  M-F  (see  Fig.  5), 
For  other  lattices,  the  procedure  is  essentially  the  same  [21, 
22].  Assume  that  an  arbitrary  point  in  the  reciprocal  space  is 
considered.  This  point  essentially  defines  a  wave  vector.  For 
the  periodicity  of  (3)  we  then  obtain  for  the  field  in  the 
original  space: 

F{r  +  R)  =  F{r)  ■  e'*'*  =  F{r)e^'^'  (5) 

where  F  stands  for  the  electric  as  well  as  for  the  magnetic 
field.  In  the  MMP  implementation  of  the  periodic  boundary 
problem  C/,  C;  are  parameters  that  characterize  the  point  in 
the  reciprocal  lattice  space.  As  a  consequence,  it  is  sufficient 
to  know  the  field  in  a  unit  cell  (as  an  equivalent 
representation  of  the  primitive  cell)  of  the  original  space.  Let 
us  call  -this  the  original  cell.  Note  that  neither  the  shape  nor 
the  location  of  the  original  cell  is  unique,  but  for  both  the 
square  and  the  hexagonal  lattice  we  may  simply  use 
quadrangular  original  cells  as  shown  in  Fig.  3  and  4. 

For  the  square  lattice,  the  relation  between  the  periodic 
constants  {Cj,  Ci)  and  the  position  in  the  IBZ  is  very  straight 
forward,  i.e.,  these  are  the  Cartesian  components  of 

the  wave  vector  k  .  For  the  hexagonal  lattice,  the  situation  is 
a  bit  more  complicated  [23-25]. 

III.  The  MMP  Solution  of  Periodic  Problems  with 
Fictitious  Boundaries 

Any  software  for  computing  band  diagrams  must  handle  both 
eigenvalue  problems  and  periodic  structures.  The  MMP 
implementation  of  MaX-1  contains  a  simple  concept  for 
handling  arbitrary  periodic  structures:  First,  the  structure  is 


Fig.  4:  The  unit  cell  of  the  photonic  crystal  with  dielectric  rods  arranged  on  a 
hexagonal  lattice. 

subdivided  into  cells  by  an  appropriate  grid  of  fictitious 
boundaries  (dashed  lines  in  Fig.  3  and  Fig.  4).  Assume  that 
the  field  in  one  of  the  infinitely  many  cells  is  known,  then, 
the  field  in  all  other  cells  is  easily  obtained  fi-om  the 
periodicity  conditions  (5),  i.e.,  the  Floquet  theorem  [7]. 

The  geometric  shape  of  the  original  cell  depends  on  the 
crystallographic  structure  (i.e.  the  crystal  symmetry),  but  it  is 
not  unique  for  a  given  crystallographic  structure  at  all, 
because  the  fictitious  boundaries  we  have  introduced,  are 
quite  ambiguous.  For  example,  in  Fig.  3  we  used  straight 
lines  between  the  circular  rods.  We  could  replace  these  lines 
by  curved,  periodic  lines  and  we  could  move  these  lines  to 
any  other  position  in  space.  Since  we  will  impose  so-called 
periodic  boundaiy  conditions  along  the  fictitious  boundaries 
of  the  original  cell,  we  have  to  minimize  the  numerical 
problems  when  we  select  the  fictitious  boundaries  in  such  a 
way  that  the  electromagnetic  field  along  them  is  as  well 
behaved  as  possible.  Therefore,  straight  lines  in  the  middle 
between  neighbor  rods  are  most  reasonable  when  the  rods  are 
circular  or  rectangular.  When  the  geometric  shape  of  the  rods 
is  more  complicated,  it  may  be  advantageous  to  use  curved 
lines. 

Once,  the  original  cell  is  isolated  by  introducing  fictitious 
boundaries,  we  can  derive  boundary  conditions  for  the  field 
along  them.  In  2D  PhCs,  the  original  cell  is  bound^  by  two 
pairs  of  parallel  lines.  For  example,  wheit  F  is  a  point  on  the 
left  border  of  the  'original  cell  in  Fig.  3,  F+e,  is  the 
corresponding  point  on  the  right  border,  where  gj 

corresponds  to  one  of  the  primitive  lattice  vectors.  Because  of 
the  periodicity,  we  obtain  from  (5): 

F(F  +  ?j)=:F(F)-e'^^'*''l^  (S') 

This  condition  holds  for  both  the  electric  and  the  magnetic 
field  in  every  point  along  the  right  boundary  of  the  original 
cell.  We  call  this  the  periodic  boundary  condition  that  is 
imposed  on  the  right  border  of  the  original  cell.  Similarly,  we 
can  introduce  a  periodic  boundary  condition  for  the  upper 
border. 
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Fig.  5:  The  band  diagram  of  Ihc  photonic  ciystal  with  dielectric  rods  on  a  square  lattice  (for  H-polarization).Thc  algorithms  used  within  the 
eigenvalue  search  procedure  arc  labeled  correspondingly. 

Having  defined  the  original  cell  and  its  periodic  boundary  we  need  at  least  two  different  expansions,  namely  a  multipole 
conditions,  one  has  to  setuptheMMP  model  of  the  scattering  expansion  (8)  or  (9)  and  Bessel  expansion  (6)  or  (7).  Note 

body  in  the  lattice  point:  We  approximate  the  field  in  each  that  the  moltipole  expansion  essentially  accounts  for  the  field 

domain  by  a  superposition  of  multipole  expansions  and  some-  scattered  at  the  inner  boundary,  whereas  the  Bessel  expansion 
times  by  additional,  analytic  solutions  of  Maxwell’s  equations  accounts  for  the  outer,  ficitious  boundaries.  This  means  that 
(in  the  frequency  domain).  The  amplitudes  or  parameters  of  the  Bessel  expansion  simulates  the  field  that  comes  from  all 
the  rcsulling  scries  expansions  are  then  computed  with  the  rods  outside  the  original  cell.  According  to  Vekua  [24],  our 
generalized  point  matching  technique,  i.c.,  b>'  minimizing  a  set  of  expansions  is  complete  in  the  sense  that  the  error  of  the 

weighted  error  fimetion  defined  along  all  natural  and  fictitious  field  is  below  an  arbitrarily  small  value  £  provided  that  the 

boundaries.  For  example,  for  the  simple  geometry  in  Fig.  4  highest  orders  are  big  enough  and  provided  that  the 
we  use  the  following  expansions:  amplitudes  (A,  B,  C,  D  in  (6)-(9) )  are  computed  correctly. 


co${n^)  +  B^J^(Kr-)sm{n(p)  (^) 

H.  =  J,(M‘)cos(n(p)  +  B^  J„(x?‘)sin(;i^)  (^) 

*t=0 

E.  =  2  cl  Hi  (w)  cos(7i^!7)  +  //^  ( w)  sin(7i  g)) 

-  — 

11=0 


where  J„  is  Bessel  function  of  order  77,  hJ  is  Hankel  function 
of  first  kind  and  order  77,  ^is  transverse  propagation  constant 
and  (r,^)  arc  polar  coordinates  wth  respect  to  the  origin  at 
the  position  of  the  corresponding  expansion.  Expansion  (6) 
(Bessel  expansion)  is  used  in  the  case  of  E-polarization  and 
expansion  (7)  is  used  in  the  case  of  H-polarization.  These 
Bessel  expansions  are  used  for  the  domain  inside  of  the 
dielectric  rod  because  these’  functions  have  no  singularity  at 
origin.  Furthermore,  these  expansions  are  sufficient  because 
the  domain  is  simply  connected.  The  background  domain  is 
not  simply  connected,  because  it  contains  a  hole.  Therefore, 


IV.  The  MMP-MAS  Eigenvalue  Solver 

For  obtaining  the  band  diagram  of  a  PhC,  it  is  necessary  to 
solve  an  eigenvalue  problem,  because  there  is  no  excitation 
like  in  scattering  problems.  This  means  that  we  only  obtain 
non-trivial  solutions  (i.e.  frequencies)  for  an  arbitrary  point  of 
the  IBZ^.e,,  for  a  given  set  of  complex  values  C/,  C2).  Thus, 
we  essentialfy  have  a  periodic  resonator  problem  to  solve.  The- 
search  of  resonance  frequencies  in  the  MMP  code  MaX-1  is 
somehow  different  from  many  other  numerical  methods 
because  MMP  uses  a  full  rectangular  system  matrix  obtained 
from  the  generalized  point  matching  technique.  For  such  type 
of  matrix  it  is  very  demanding  to  obtain  accurate  results  while 
avoiding  problems  with  the  condition  number  [25].  Note  that 
condition  number  problems  are  especially  crucial  when  one  is 
solving  eigenvalue  problems.  If  this  is  not  properly  done,  one 
can  obtain  a  "noisy"  behavior  near  the  eigenvalues  and  this 
can  heavily  disturb  the  numerical  eigenvalue  search 
procedures.  However,  the  standard  MMP  eigenvalue  search 
procedure  first  defines  a  real  valued,  positive  eigenvalue 
search  function 
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A\e) 


(10) 


where  e  is  the  eigenvalue  (i.e.  the  resonance  frequency),  E  is 
the  weighted  residual,  and  A  is  an  amplitude  that  may  be 
retrieved  from  any  field  component  in  a  specific  test  point  (or 
an  integral  over  some  field  profile).  For  the  band  gap 
computation,  it  is  most  reasonable  to  define  A^  as  the  total 
electromagnetic  energy  within  the  original  cell.  According  to 
(10)  the  desired  eigenvalues  are  characterized  by  the  minima 
of  the  search  function  77.  Analyzing  the  shape  of  77  near  the 

minima  provides  additional  information  on  the  accuracy  of 
the  solution. 

Although  more  reliable  results  are  obtained  when  the 


T  -  - —  it  K-vccior  in  me  i 

point)  using  the  fictitious  excitation  in  a  random  and  symmetric  position 
respectively. 


amplitude  is  defined  by  an  appropriate  integral,  the  definition 
in  one  or  a  few  test  points  is  sufficient  for  most  cases.  Since 
the  numerical  intergration  may  be  time-consmning,  one 
usually  prefers  the  simpler  test  point  method  However,  it  is 
important  to  note  that  the  definition  of  the  search  function  is 
not  unique.  By  defining  different  search  functions,  one  can 
gain  even  more  intrinsic  information  providing  a  good  error 
estimation  and  for  validation  puiposes.  As  depicted  in  Fig.  6, 
even  for  a  single  model  (fixed  amplitude  definition  and  fixed 
multipole  expansion),  one  can  address  the  different  minima  of 
the  same  eigenvalue  search  functions  simply  by  rearranging 
the  columns  of  the  MMP  matrix.  In  feet,  in  the  Givens  update 
algorithm  [25],  which  was  used  for  solving  the  MMP  matrix 
equation,  the  last  expansion  somehow  plays  the  role  of  an 
excitation.  When  it  happens  that  the  spatial  symmetry  of  such 
excitation  is  not  contained  in  the  symmetry  of  the  searched 
eigenmodc,  this  mode  will  not  be  "excited”,  hence,  the 
corresponding  minimum  of  the  eigen-value  search  function  is 
suppressed.  Although,  it  may  be  desirable  to  suppress  some 
modes  in  ^plications  where  not  all  modes  must  be 
considered,  this  is  usually  inconvenient  for  the  automatic 
computation  of  the  complete  band  structure.  We  therefore 
look  for  an  alternative  technique. 

Remember  that  we  have  introduced  fictitious  boundaries  for 
handling  the  periodic  problem.  Similarly,  we  now  can 
introduce  a  fictitious  excitation  that  is  defined  in  such  a  way 
that  all  modes  are  excited  (Fig.  7).  This  concept  mimics  the 
measurement  of  resonance  frequencies,  where  one  always 
needs  an  excitation  of  the  resonator  and  a  test  point  (or  port) 
where  the  signal  is  measured.  By  sweeping  the  frequency  of 
the  exitation,  the  peaks  of  the  amplitude  A  in  the  test  point 
CM  be  readily  assigned  to  the  resonance  frequencies  of  the 
different  modes.  This  procedure  was  first  introduced  by  the 
Method  of  Auxiliary  Sources  (MAS)  [26]  and  a  similar 
method  was  used  by  Sakoda  [27].  Finally,  the  method  was 
adapted  to  MMP  by  Moreno  [28].  MAS  uses  eigenvalue 


fictitious  excitation. 
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search  functions  such  as  the  energy  density  at  the  test 
point  are  used.  The  eigenvalues  are  then  obtained  from  the 
maxima  of  p.  The  analysis  of  \i  near  the  maxima  has  yielded 
a  strange  "double  peak"  phenomenon  that  disturbs  the 
numerical  search  procedure.  The  standard  MMP-MAS 
eigenvalue  solver  searches  for  minima  of  the  eigenvalue 
search  function  i*e.,  one  obtains  "twin 

minima"  instead  of  double  peaks,  as  shown  in  Fig.  8.  The 
"double  peak"  phenomenon  and  the  "twin  minima"  are  caused 
by  a  very  sharp  peak  of  the  residual  E  at  the  correct 
eigenvalue  position.  Note  that  this  peak  is  not  obtained  in  the 
standard  MMP  approach  without  fictitous  excitations.  Of 
course,  the  residual  peak  may  also  be  used  for  defining  the 
eigenvalues.  Since  these  peaks  are  extremely  sharp,  it  is  very 
likely  that  one  of  the  eigenvalues  is  missed  by  the  rough 
search  routine  that  searches  for  all  eigenvalues.  In  order  to 
overcome  these  problems,  one  can  define  more  complicated 
eigenvalue  search  functions  t]  as  proposed  in  Fig.  8.  This  allows 
one  to  suppress  the  double  peak  phenomenon.  Unfortunately, 
one  may  encounter  numerical  underflow  problems  in  some 
applications.  Therefore,  the  current  MaX-l  eigenvalue  solvers 


uses  three  different  "competing"  eigenvalue  search  functions: 
1)  A  complicated  one  vwth  user-definable  exponent  n,  2)  the 
inverse  of  the  amplitude,  and  3)  the  proper  residual.  Using  all 
of  these  three  functions,  the  code  is  capable  to  detect  the 
correct  locations  of  the  eigenvalues.  An  alternative  to 
overcome  the  twin  minima  problems  is  the  introduction  of 
"fictitious  losses"  that  smoothen  the  resulting  search  function 

?7- 

Since  one  often  considers  a  broad  frequency  range,  it  is  not 
reasonable  to  find  the  eigenvalues  by  plotting  the  eigenvalue 
search  function  over  the  entire  range  with  a  very  high 
resolution.  It  is  much  more  efficient  to  subdivide  the  search 
process  into  two  steps:  1 )  Rough  detection  of  all  eigenvalues 
and  2)  fine  search,  i.e.,  accurate  computation  of  the 
eigenvalues.  The  first  step  seems  to  be  trivial  as  soon  as  the 
problems  mentioned  above  have  been  solved.  The  second  step 
requires  a  fast  minimum  search  procedure  for  real  flincticxis. 
The  algorithm  used  In  MaX-1  is  mainly  based  on  a  parabolic 
interpolation  because  the  search  function  near  the  minima  is 
usually  almost  parabolic  provided  that  the  double  peak 
phenomenon  has  been  removed. 


Fig.  9:  The  algorithm  for  the  band  structure  computation  using  MMP. 


!  from  gamma  to  sc,  left  part 

read  win  sc|uare_gjc_m_gjTOcJdle.vvin 

draw'windcMrl  Xm 

read  turction  m. fun 

drowfunctlonl  2  1 

read  function  xJun 

draw  function  121 

rirowtext  173  400 1  "x" 

drewvtext  349  400 1  "m" 

ciraw  text  0  400 1  “gamma** 

drawtod  560  400 1  "gamma" 

write  function  *  header  0  2 

write  function  /  teact  **11* 

wrrte  function  /  teact  “bratouin  zone" 

write  function  /  text  **18*:“ 


write  fun  ! 

read  mmp  bandOOl  .mmp 
set  eigenvalue  fine  1 
end  pre  preressing  cfirectives 
! 

-►  set  period  constant  1 57079B.3  0  0  0 
set  var  41 

write  function  header  0  2 
wrltefunctlon  /  tod  "rf 
writefunction  /  text  "fcrinoLiin  zone" 
write  function  /  tod  Tafc" 
end  outer  loop  pre  drectwes 


mmp  solve 
writefun  /  var 
writefun/lrerea 
increase  period  constert 
increase  ver  -1 
end  inner  loop  directVes 


39269.9075  0 


writefun  ! 
read  function/ 

process  fu  ndion  clivC^/2.3c1 4) 
drawfunctioni  3  1 
Increase  eigenvalue  1 
end  outer  loop  post  directives 


->i 


read  dir  sc|uare_5jf_m_gjTiiddle_2.dir 


00 


Fig.  10:  The  algorithm  for  band  diagram  computation  written  in  the  MaX-l 
script  language. 
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Table  I 

Convergence  Characterjstics,  Computation  for  l  ”  and  6"*  Eigenfrequency  at  X  Point  of  IBZ 


Number  of 
unknowns 

Eigenfrequency  1 

Eigenfrequency  6 

Frequency  (Hz) 

Error  (%) 

Field  mismatch.  (%) 

Frequency  (Hz) 

Error  (%) 

Field  mismatch.  (%) 

20 

1.0223S85el4 

1.473 

9.620748e-0 

2.3465308el4 

4.194 

2.306028C+1 

36 

I.0095955el4 

0.206 

4.204639e-0 

2.2567438el4 

0.207 

4.449937C+0 

52 

1.0078338el4 

0.032 

0.2881 15e-0 

2.2S12072el4 

0.039 

0.824365e43 

94 

1.0074678cl4 

0.005 

4.465747C-2 

2.251942IeI4 

0.006 

0.128581e-0 

164 

1.0075153cl4 

0.000 

4.72972  Ie-7 

2.2520785cl4 

0.000 

3.551224C-6 

Having  a  closer  look  to  topical  band  diagrams  (Fig.  3),  we  see 
different  situations  which  can  cause  problems  for  both  the 
rough  search  and  the  line  search.  Mainly  at  the  T  and  the  M 
point  we  usually  observe  degenerate  modes.  Furthermore,  we 
have  areas  with  almost  degenerate  modes  and  points  i^ere 
different  lines  seem  to  cross  each  other,  where  the  modes  are 
(accidentally)  degenerated.  When  the  rough  search  is 
performed  to  degenerate  points,  it  usually  cannot  detect  all 
modes  involved.  Even  if  the  search  procedure  is  started  in  a 
close  vicinity  to  such  degeneracies,  it  will  be  too  time- 
consuming  to  iterate  into  all  eigenmodes.  In  order  to 
overcome  these  problems,  it  is  reasonable  to  start  a  rough 
search  in  a  domain  where  all  eigenvalues  could  be  easily 
tracked  down  (e.g.  the  interval  between  T  and  X  in  the  band 
diagram  of  Fig.  5).  Once  this  has  been  done,  one  can  trace 
each  eigenvalue  by  moving  a  small  step  either  to  the  left  or 
right  side  within  the  band  diagram,  and  repeating  this 


Fig.  1 1 :  The  band  diagram  of  the  photonic  crystal  with  dielectric  rods  and  square 
lattice,  H-polartzadon,  the  first  6  modes. 


Fig.  12:  The  band  diagram  of  the  photonic  crystal  vdth  dielectric  rods  and 
square  lattice,  E-polarization,  the  first  6  modes. 


procedure  until  the  border  of  the  diagram  is  reached.  For  each 
such  step,  only  a  ftne  search  must  be  performed.  Depending 
on  1)  the  desired  accuracy,  2)  the  step  size,  and  3)  special 
properties  of  the  model,  several  iterations  are  required.  The 
number  of  iterations  could  be  drastically  reduced  when  using 
the  Eigenvalue  Estimation  Technique  ^ET)  implemented  in 
MaX-1  [19].  This  technique  uses  the  information  of  previous 
eigenvalue  solutions  for  the  extrapolation  of  the  current 
eigenvalue’s  search  interval.  Typically,  4-12  iterations  per 
step  are  sufficient  for  obtaining  an  eigenvalue  with  a  high 
precision.  For  example,  for  tracing  the  first  mode  in  Fig.  3, 
280  search  steps  were  performed  and  5  iterations  per  step 
were  required  in  the  average. 

V.  Automatic  Eigenvalue  Search 

Referring  to  e.g.  Fig.  5  a  standard  band  diagram  consists  of 
three  different  intervals  corresponding  to  the  three  sides  of 


Fig.  13:  The  band  diagram  of  the  photonic  crystal  with  dielectric  rods  and 
hexagonal  lattice,  H-polarization,  the  first  6  modes. 


Fig.  14:  The  band  diagram  of  the  photonic  crystal  with  dielectric  rods  and 
hexagonal  lattice.  E-polarization,  the  first  6  inodes. 
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the  IBZ.  When  the  rough  search  is  started  somewhere  in  the 
middle  of  such  an  interval  (e.g.  in  the  area  between  F  and  X 
in  the  band  diagram),  it  must  be  repeated  three  times.  After 
each  rough  search  the  fine  search  must  be  repeated  for  each 
obtained  eigenvalue  and,  finally,  the  fine  search  routine  must 
run  for  each  eigenvalue  once  towards  the  left  and  once 
towards  the  right  side  of  the  band  diagram,  as  depicted  in  Fig. 
5.  MaX-1  contains  a  script  language  that  allows  one  to  define 
complicated  procedures  such  as  the  search  procedure  mentioned 
above.  The  set  of  MaX-1  directives  for  the  automatic 
generation  of  a  band  diagram  from  the  point  in  the  middle 
between  F  and  X  to  the  F  point,  is  given  in  Fig.  10,  and  the 
complete  algorithm  for  this  procedure  is  given  in  Fig.  9.  It  is 
obvious  that  the  algorithm  is  not  simple  and  the  overall 
procedure  relics  on  fast  computer  resources. 

VI.  Numerical  Verification 

We  have  applied  MMP  to  various  PhC  lattices.  Internal  tests 
show  excellent  convergence.  Therefore  high  accuracy  may 
easily  be  obtained.  Table  I  shows  the  MMP  estimate  of  the 
mismatching  errors  along  the  boundary  for  the  model 
outlined  in  Fig.  3  with  different  maximum  orders  of  the 
multipolcs  and  Bessel  expansions,  i.e.,  with  different  numbers 
of  unknowns.  Note  that  the  computation  time  typically  is 
proportional  to  the  cube  of  the  number  of  unknowns  because 
wc  use  a  brute-force  fiill  matrix  solver  (Givens  update 
scheme).  Despite  of  this,  *  the  computation  time  remains 
reasonably  short  because  the  matrices  obviously  are  much 
smaller  than  the  matrices  used  in  other  methods.  For  example 
Fig.  11  was  obtained  with  3  rough-search  routines,  100 
frequency  slcp.s  each.  The  total  number  of  1656  plotted  points 
required  were  then  computed  with  8280  MMP  evaluations  of 
i.e.  approximately  5  iterations  per  point  in  the  diagram 
were  performed.  The  total  calculation  time  was  40  minutes  on 
a  Pentium  4,  2GHz.  Because  of  the  excellent  convergence,  we 
also  can  estimate  the  accuracy  of  the  eigenvalues  by 
comparing  them  with  a  very  accurate  MMP  model.  As  one 
can  sec  from  Table  I,  one  only  obtains  one  more  digit  when 
doubling  the  number  of  unknowns. 

In  order  to  validate  this  algorithm,  several  calculations  were 
performed  and  results  were  compared  with  the  results  of  MPB 
package  developed  at  the  MIT  [29].  For  the  PhC  with  square 
lattice  and  dielectric  rods  (Fig.  3),  a  band  (fiagram  calculation 
was  performed  for  different  field  polarizations  and  the  results 
arc  given  in  Fig.  11  (H-polarization)  and  Fig.  12  (E- 
polarization).  The  results  for  the  hexagonal  lattice  case  (Fig. 
4),  arc  depicted  in  Fig.  13  (H-polarization)  and  Fig.  14  (E- 
polarization).  These  two  types  of  PhC  rely  on  the  same  lattice 
data:  A  dielectric  rod  with  radius  r  =  0.3a  and  a  dielectric 
constant  of  £  =  11.56,  the  lattice  is  embedded  in  air  and  the 
lattice  constant  is  a  =  lO^(m)*  From  Figs.  1 1-14  we  deduce  a 
perfect  agreement  with  the  MPB  results  documented  in  [29]. 

VII.  Conclusion 

Wc  have  presented  an  efficient  method  for  band  structure 
calculation  for  2D  dielectric  PhCs.  In  this  framework  a  fully 


automatic  algorithm  was  developed  and  evaluated  along 
several  examples.  The  eigenvalue  searching  procedure  in  the 
frequency  domain  has  been  performed  using  a  fictitious 
excitation.  Optimal  eigenvalue  search  functions  have  been 
found  while  evaluating  the  total  eigenvalue  spectrum  for  k- 
valucs  at  three  preferable  points  on  the  IBZ.  The  three 
resulting  sets  of  eigenvalues  are  evolved  into  a  foil  band 
diagram  using  a  highly  efficient  Eigenvalue  Estimation 
Technique  (EET).  The  overall  algorithm  performs  photonic 
band  diagram  calculations  at  a  very  high  level  of  accuracy 
and  at  reasonable  computational  costs.  This  algorithm  is 
easily  extendable  for  applications  involving  localized  defect 
mode  analysis  [30],  various  PhC  defect  waveguide  types 
(supcrcell  approach  [31])  and  photonic  waveguide 
discontinuities  [31],  as  well. 
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Abstract 

Delay  lines  are  convenient  circuit  elements  used 
to  introduce  delay  between  circuit  board  components 
to  achieve  required  timing,  Seipentine  or  meander 
lines  arc  the  most  common  of  delay  lines.  These  lines 
introduce  delay  but  also  introduce  spurious 
dispersion  that  makes  the  signal  appear  as  if  it  is 
arriving  earlier  than  expected.  The  cause  of  such 
spurious  speed-up  or  skew  is  analyzed  qualitatively. 
Previous  work  found  that  owing  to  the  periodicity 
inherent  in  the  serpentine  line  structure,  the  crosstalk 
noise  accumulates  synchronously,  thus  creating  a 
higher  potential  for  triggering  false  logic.  Numerical 
simulations  arc  performed  using  the  Finite-Difference 
Time-Domain  (FDTD)  method  to  corroborate  the 
qualitative  prediction  with  physical  behavior.  Based 
on  the  understanding  of  the  coupling  mechanism  in 
periodic  serpentine  lines,  a  qualitative  prediction  can 
be  made  of  the  behavior  of  novel  delay  lines  such  as 
the  spiral  line.  A  new  delay  line,  the  concentric  Cs 
delay  lines,  is  introduced.  The  design  of  the  new  line 
is  based  on  forcing  the  crosstalk  noise  to  spread  over 
time,  or  to  accumulate  asynchronously,  thus 
enhancing  the  integrity  of  the  received  signal. 

Key  Words:  Delay  Lines,  Serpentine  Lincs,"Meander 
Lines,  Spiral  Lines,  FDTD,  Numerical  Simulation, 
Crosstalk  Noise. 


I.  Introduction 

It  can  be  argued  that  timing  problems  are 
amongst  the  most  serious  plaguing  high-speed  digital 
circuit-board  design.  As  the  clock  speed  increases, 
the  wavelength  shrinks.  For  instance,  for  a  clock 
speed  of  1.5  GHz,  the  pulse  harmonics  containing 
sufficient  energy  reach  beyond  10  GHz,  making  a 
typical  motherboard  or  daughter  cards  electrically 


large.  As  a  consequence,  a  delay  line  used  to 
introduce  precise  delay  (based  on  the  length  of  the 
line  and  board  material)  between  circuit  board 
elements  can  no  longer  be  considered  to  have  TEM 
or  quasi-TEM  propagation  behavior,  and  to  be 
electromagnetically  isolated  from  neighboring  lines 
that  fall  within  its  proximity.  With  the  timing  budget 
shrinking  to  the  picosecond  regime,  any  non- 
predictable  behavior  of  delay  lines  can  potentially 
cause  a  timing  imbalance  at  the  receiver  end. 

Two  mechanisms  are  employed  to  achieve 
required  signal  delay  between  circuit  components. 
The  first  is  achieved  through  internal  electronic 
circuitry.  The  second  mechanism,  which  is  the  most 
common  and  least  expensive,  is  achieved  through 
meandering  a  transmission  line  as  shown  in  Fig.  1. 
The  meandered  line,  commonly  referred  to  as  the 
serpentine  line,  consists  of  a  number  of  closely 
packed  transmission  line  segments.  The  only 
objective  from  meandering  the  line  is  to  achieve  high 
density  (of  transmission  line)  per  square  inch  of 
circuit  board  while  obtaining  signal  delay  that  is 
directly  proportional  to  the  length  of  the  line. 

When  serpentine  lines  arc*  used  in  high-speed 
digital  circuit  applications,  the  time  delay  through  a 
single  serpentine  line  can  be  much  longer  than  the 
rise  time  of  the  signal  pulse.  Under  such  conditions, 
serpentine  lines  have  been  found  to  introduce  a 
dispersion  that  makes  the  signal  appear  as  if  it  is 
arriving  earlier  than  would  be  expected  based  on  the 
exact  electrical  length  of  the  line  [1].  This  spurious 
speed-up  in  the  signal  can  be  characterized  as  a  type 
of  dispersion.  This  dispersion  is  primarily  caused  by 
the  topology  of  the  line  and  is  directly  related  to  the 
crosstalk  between  the  adjacent  transmission  line 
sections.  More  specifically,  the  dispersion  is  related 
to  two  fundamental  geometrical  parameters:  the  first 
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Fig.  1.  A  matched  seven-section  serpentine  delay 
line.  The  Characteristic  impedance  of  each  identical 
section  is 

is  the  length  of  each  serpentine  section,  and  the 
second  is  the  spacing  between  adjacent  sections, 
which  influences  the  capacitive  and  inductive 
coupling  [1],  [2]. 

Based  on  these  findings,  Wu  and  Chao 
introduced  the  flat  spiral  delay  line  to  force  the 
crosstalk  noise  to  accumulate  asynchronously  [3]. 
The  flat  spiral  line  is  distinguished  fiom  the  classical 
spiral  delay  line,  which  requires  three-dimensional 
topology  and  therefore  cannot  be  printed  on  a  single 
board  layer.  In  this  work,  the  flat  spiral  line  will  be 
referred  to  as  the  spiral  line. 

The  spiral  delay  line  is  considered  a  significant 
improvement  over  the  serpentine  line.  In  [3],  the 
motivation  behind  the  spiral  line  was  introduced  and 
analysis  was  given  based  on  the  wave-tracing 
technique  which  did  not  take  into  account  higher- 
order  mode  coupling  and  comer  effects.  Also  in  [3], 
quantitative  analysis  was  also  performed  to  include 
the  effects  of  multiple  and  feedback  couplings 
between  lines.  To  incorporate  high-order  effects  into 
the  modeling  of  serpentine  line,  especially  the 
coupling  between  the  non-orthogonal  lines,  the  three- 
dimensional  Finite-Difference  Time-Domain  (FDTD) 
method  was  used  to  analyze  serpentine  and  spiral 
lines  [4].  In  a  subsequent  work,  the  FDTD  method 
was  used  to  analyze  the  spiral  line  for  Gaussian  pulse 
excitation  [5].  More  recently,  the  Method  of 
Moments  (MoM)  was  used  to  present  a  full-wave 
model  for  the  serpentine  line  including  physical 
losses  [6]. 

In  this  work,  we  use  the  same  qualitative  wave¬ 
tracing  methodology  that  was  adopted  in  [2]  to 
analyze  the  spiral  line.  We  show  that  the  wave¬ 
tracing  methodology  can  only  give  qualitative 
analysis  of  the  line  and  that  a  full-wave  model  is 
needed  to  accommodate  the  effect  of  comers  and 


coupling  between  non-orthogonal  segments.  We  use 
the  FDTD  method  to  analyze  the  spiral  line  and  show 
that  more  inclusive  performance  can  be  predicted. 
Based  on  the  qualitative  insight  gained  from  the 
analysis  of  the  serpentine  and  spiral  line,  we 
introduce  a  new  delay  line  that  we  refer  to  as  the 
concentric  Cs  delay  line.  Qualitative  analysis  of  the 
new  lines  is  given  followed  by  full-wave  three- 
dimensional  simulation. 

The  organization  of  this  paper  is  as  follows.  In 
section  2,  the  mechanism  of  coupling  between  two 
tran^ission  lines  is  introduced  with  application  to 
seipentine  lines.  Sections  3  and  4  present  a 
qualitative  analysis  followed  by  full-wave  numerical 
simulation  of  practical  real-world  structures  for  the 
serpentine  and  spiral  lines,  respectively.  Section  5 
introduces  the  Concentric  Cs  delay  lines.  We  note 
that  the  FDTD  method  will  be  used  purely  as  a 
numerical  simulation  tool  without  giving  details  of 
the  simulation  parameters  involved  such  as  the  cell 
size  and  time  step,  ...etc.  The  paper  concludes  with 
critical  observations  that  can  be  used  in  the  design  of 
predictable-delay  lines. 

n.  Weak  Coupling  Crosstalk 

a.  Analysis 

A  typical  serpentine  delay  line  is  composed  of 
transmission  line  sections  closely  packed  as  shown  in 
Fig.  1.  Let  us  isolate  two  adjacent  sections  as  shown 
in  Fig.  2,  After  examination,  we  notice  that  the  two 
isolated  sections  resemble  the  simplest  case  of  two 
parallel  transmission  lines  [2].  After  ignoring  higher- 
order  effects,  the  two-segment  serpentine  line  shown 
in  Fig,2  is  equivalent  to  the  two  transmission  lines 
matched  at  both  ends  and  shown  in  Fig.  3.  Once  this 
observation  is  made,  using  conventional  transmission 
line  analysis,  one  can  predict  the  crosstalk  induced  on 
the  second  line  due  to  the  launched  signal  on  the  first 
line. 


n 


•  •  • 


Fig.  2.  Segment  of  a  serpentine  line. 

Considering  Fig.  3,  the  near-end  crosstalk, 
shown  as  voltage  V2,  is  proportional  to  the  mutual 
capacitance  and  mutual  inductance  characterizing  the 
two  lines.  The  coupling  coefficient  is  given  by 
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(1) 


“  Ac,, 

where  the  crosstalk  at  the  near  end  is  given  by 


To  demonstrate  these  two  scenarios,  we  use  the 
FDTD  method  to  simulate  the  behavior  of  two 
parallel  strip  lines  (similar  topology  to  Fig.  3)  with 
cross  section  shown  in  Fig.  4. 


where  Cn  and  Ln  are  the  self  capacitance  and 
self  inductance,  respectively,  of  the  lines,  and  Cn 
and  L]2  arc  the  mutual  capacitance  and  mutual 
inductance,  respectively,  of  the  coupled  lines  [I].  Let 
us  assume  that  the  rise  time  is  smaller  than  the  line 
delay,  tj.  If  a  step  function  is  launched  on  line  1,  then 
the  crosstalk  at  the  near  end  of  line  2  will  be  a  pulse 
with  a  duration  twice  the  line  delay  (i.e.,  equivalent 
to  the  round-trip  time  tj).  The  crosstalk  at  the  far  end 
has  a  much  smaller  duration  and  it  equals  zero  when 
the  capacitive  and  inductive  coupling  coefficients  are 
equal  (as  in  the  case  when  the  medium  is 
homogeneous  as  in  strip  lines)  [l]-[2].  Since  the 
pulse  on  the  active  line  propagates  to  the  right,  the 
voltage  at  the  near  end  of  the  passive  line  is  in  effect 
due  to  a  leftward  propagating  wave.  Similarly,  the 
voltage  at  the  far  end  of  the  passive  line  is  the 
accumulation  of  a  rightward  propagating  wave 
(which  is  zero  in  the  case  of  homogeneous  medium. 


Fig.  3.  Parallel  transmission  lines  with  matched 
terminations.  A  pulse  transmitted  on  the  upper  line 
generates  crosstalk  on  the  lower  line  that  propagates 
in  a  manner  equivalent  to  the  case  of  the  serpentine 
line  segment  shown  in  Fig.  2. 

If  line  1,  in  Fig.  3,  is  excited  by  a  pulse  starting 
at  t=  to  and  of  duration  T  assumed  to  be  shorter  than 
the  length  of  the  line  td,  then  the  near-end  crosstalk 
will  consist  of  two  opposite  polarity  pulses  separated 
by  2td-T,and  each  of  duration  T.  This  can  be  shown 
using  the  principle  of  superposition.  By 
decomposing  the  finite  duration  input  pulse  into  two 
opposite  polarity  step  functions  separated  by  T,  the 
near-end  crosstalk  will  be  the  sum  of  two  pulses.  The 
first  pulse  will  have  positive  polarity  and  is  excited  at 
time  to.  The  second  pulse  will  have  negative  polarity 
and  will  start  at  time  to+T.  The  sum  of  the  two  pulses 
yields  a  waveform  consisting  of  two  pulses,  each  of 
duration  T  and  separated  by  a  time  interval  of  width 
2td-T. 


W=0.}mm,  L-Jmm,  H=0,4mm,  5 


Fig.  4.  Strip  line  cross  section. 

Figure  5  shows  the  near-end  voltages  on  the 
active  and  passive  lines  when  the  excitation  is  a  step 
function.  The  far-end  voltage  at  the  active  line  is  also 
shown  in  Fig.  5,  One  can  observe  that  the  duration  of 
the  near-end  crosstalk  is  T. 

Figure  6(a)  shows  the  near-end  voltage  on  the 
active  and  passive  lines  when  the  excitation  is  a  pulse 
of  duration  250  picoseconds.  We  observe  that  the 
near-end  passive  line  experiences  two  opposite 
polarity  pulses,  each  of  the  same  duration  as  the 
original  pulse.  This  performance  is  in  full  agreement 
with  our  earlier  prediction.  For  completion,  we  show 
the  far-end  voltage  on  the  active  and  passive  lines  in 
Fig.  6(b).  Notice  that  the  bump  appearing  in  Fig.  5  is 
due  primarily  to  the  numerical  dispersion  of  the 
FDTD  method  and  the  physical  dispersion  caused  by 
the  propagation  in  the  strip  line  configuration. 

When  forming  a  serpentine  line  by  closely 
connecting  identical  transmission  line  sections,  a 
similar  crosstalk  mechanism  to  that  discussed  above 
takes  effect,  except  now  one  should  also  account  for 
coupling  between  non-orthogonal  lines  and  the 
coupling  between  orthogonal  line  segments.  Let  us 
consider  the  serpentine  line  composed  of  seven 
transmission  line  segments  (shown  in  Fig.  1). 
Suppose"  Vjn  is  a  step  voltage  with  a  magnitude  of  one 
volt.  - 

Once  this  voltage  waveform  is  launched  at  the 
transmitter  end,  a  voltage  of  magnitude  1/2  appears 
on  line  1.  The  voltage  at  line  1  will  induce  leftward 
traveling  waves  at  the  near  ends  of  lines  2,  4  and  6* 
These  signals  continue  to  travel  on  lines  3,5  and  7 
respectively.  Notice  that  the  launched  signal  does  not 
excite  any  forward  propagating  waves  on  any  of  the 
lines  as  per  the  discussion  given  above.  (Throughout 
this  work,  forward  propagating  cross  talk  refers  to 
crosstalk  propagating  towards  the  receiver.  Similarly, 


Ramahi:  Analysis  of  Conventional  and  Novel  Delay  Lines 


184 


“backward  propagating  crosstalk  refers  to  crosstalk 
propagating  towards  the  source.) 


0.0  OJ  a4  0.S  M  1JI  t.3 

time  (nsnoMcontQ 

Fig.  5.  Wave  forms  on  the  active  and  passive  lines 
due  to  a  step  function  excitation.  VI  is  the  input  pulse 
at  the  active  line,  V2  is  the  near-end  voltage  at  the 
passive  line,  and  V3  is  the  far-end  voltage  at  the 
active  line. 

Assuming  that  multiple  coupling  is  negligible 
(multiple  coupling  refers  to  the  feedback  crosstalk 
between  lines),  equation  (1)  can  be  used  to  find  the 
magnitudes  of  the  crosstalk  induced  on  any  line  by 
the  main  signal  as  it  propagates  on  any  other  line. 
The  coupling  between  any  two  lines  m  and  n,  where 
m  n,  can  be  written  as 


where  Cnun  and  Lmm  are  the  self  capacitance  and 
self  inductance,  respectively,  of  line  m,  and  Cmn  and 
Lnn  are  the  mutual  capacitance  and  mutual 
inductance,  respectively,  of  the  coupled  lines  m  and 
n.  Let  us  assume  that  the  rise  time  is  much  smaller 
than  the  pulse  width.  Let  us  consider  the  leftward 
traveling  wave  excited  at  line  2,  which  has  a 
magnitude  of  ki/2,  and  duration  2td.  Furthermore,  this 
wave  arrives  at  the  receiver  end  ahead  of  the  main 
signal  by  2td.  Similarly,  backward  traveling  waves 
will  appear  on  lines  4  and  6  having  duration  of  21^ 
and  magnitudes  k3/2  and  ks/2  respectively.  These 
signals  will  arrive  ahead  of  the  main  signal  by  4td  and 
6td,  respectively. 

When  the  niain  signal  starts  to  propagate  along 
line  2,  it  induces  forward  propagating  crosstalk  on 
lines  3,  5  and  7.  The  magnitudes  of  these  crosstalk 
signals  are  ki/2,  k^/l  and  respectively.  These 


signals  arrive  ahead  of  the  main  signal  by  2td,  and 
6td,  respectively.  It  is  important  to  realize  that  as  the 
main  signal  hops  from  one  line  to  the  next,  the 
crosstalk  induced  is  additive.  In  other  words,  the 
crosstalk  noise  accumulates  synchronously.  This 
process  of  coupling  continues  until  the  main  signal 
arrives  at  the  receiver.  In  fact,  by  adding  the  effect  of 
all  induced  crosstalk,  the  waveform  obtained  at  the 
receiver  end  becomes  a  laddering  waveform,  where 
each  of  the  ladder  levels  is  of  duration  2td. 


time  (nanosecond) 

(a) 


Fig,  6.  Wave  fonrw  on  the  active  and  passwe  sections 
of  a  two-line  system  obtained  using  the  FDTD 
method,  (a)  Near-end  voltage  (NE-a  =  Near  end  at 
active  line,  NE-p  =  Near  end  at  passive  line),  (b)  Far- 
end  voltage  (FE-a  =  Far  end  at  active  line,  FE-p  =  Far 
end  at  passive  line). 

Therefore,  when  the  seven-section  serpentine 
line  (Fig.  1)  is  excited  widi  a  step  function  starting  at 
to,  we  can  expect  the  receiving  end  waveform  to  have 
a  laddering  waveform.  When  the  excitation  is  a  finite 
duration  pulse,  then  the  shape  of  the  received  signal 
will  depend  on  whether  T  is  shorter  or  longer  than  the 
round  trip  2td.  When  T  is  longer  than  2td,  the 
qualitative  behavior  of  the  received  signal  is  depicted 
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in  Fig,  7(a),  showing  degradation  of  the  pulse 
magnitude  once  the  main  signal  is  received.  On  the 
other  hand,  when  T  is  shorter  than  2td,  the  received 
waveform  consists  of  a  train  of  pulses,  each  of 
duration  T  and  separated  by  2td  -  T,  as  shown  in  Fig. 
7(b). 

TIic  wave  tracing  analysis  adopted  here  can  be 
further  generalized  to  serpentine  lines  composed  of 
many  sections.  Tt  can  be  shown  that  for  a  general 
serpentine  delay  line  consisting  of  2N+1  sections,  a 
laddering  waveform  is  generated  and  it  arrives  before 
the  main  signal.  The  laddering  wave  includes  N 
ladders;  each  ladder  is  of  width  2td  and  level 
for  i=  I  to  N. 

In  light  of  the  above  analysis,  several 
observations  can  be  made; 

1.  The  laddering  wave  creates  the  impression 
that  the  pulse  is  arriving  earlier  than 
intended.  In  reality,  however,  the  waveform 
arriving  at  the  receiver  is  composed  of  two 
signals:  the  crosstalk  laddering  waveform 
(noise)  and  the  unadulterated  signal. 

2.  The  crosstalk  accumulates  synchronously  at 
the  receiving  end.  The  synchronous 
accumulation  is  due  to  the  periodicity  of  the 
serpentine  line  structure. 

3.  The  highest  level  of  the  laddering  waveform 
can  reach  several  times  the  level  of  crosstalk 
induced  between  two  adjacent  lines.  In  fact, 
the  cumulative  magnitude  can  be  higher  than 
the  level  of  the  main  signal. 

4.  The  magnitude  of  each  ladder  is 
proportional  to  the  capacitive  and  inductive 
coupling  between  the  lines. 

5.  The  spreading  of  the  laddering  wave,  i.e., 
the  width  of  each  ladder,  is  directly 
proportional  to  the  length  of  the  serpentine 
line  segments. 

6.  When  the  pulse  duration  is  greater  than  2td, 
where  t^,  is  the  delay  through  Tone  section, 
the  cumulative  effect  is  less  pronotmted, 
producing  ladders  of  short  duration.  The  net 
effect  is  that  the  pulse  arrives  at  the 
receiving  end  as  a  dispersed  wave  with  an 
apparent  rise  time  longer  than  that  of  the 
original  signal.  In  such  case,  the  serpentine 
line  can  be  interpreted  as  a  low-pass  filter. 

7.  When  the  pulse  duration  is  smaller  than  2td, 
where  td,  is  the  delay  through  one  section, 
the  received  waveform  consists  of  a  train  of 
pulses,  each  of  duration  T  and  separated  by 
2ta-T. 


b.  Numerical  Experiments 

Several  assumptions  were  made  in  the  wave 
tracing  analysis  presented  above.  These  were: 

1 .  Negligible  non-adjacent  line  coupling. 

2.  The  backward  propagating  crosstalk  was 
assumed  to  be  zero.  In  the  general  case  of 
non-homogeneous  media,  there  is  some 
backward  induced  crosstalk. 

3.  The  small  transmission  line  segments  that 
connect  the  longer  sections  (orthogonal 
sections)  have  been  assumed  to  have  zero 
delay  (zero  physical  length). 

4.  Negligible  multi-modal  propagation. 

Considering  the  above  simplifications,  the 
qualitative  model  discussed  earlier  serves  only  to 
give  an  understanding  of  the  primary  crosstalk 
contributors.  For  a  more  accurate  prediction  of  the 
performance  of  the  serpentine  line,  the  three- 
dimensional  FDTD  method  will  be  used  since  it  fully 
integrates  the  constraints  listed  above.  ^ 

The  first  design  tested  is  that  of  the  seven- 
section  serpentine  line  disaisscd  earlier  (see  Fig.  1). 
A  cross  section  of  two  adjacent  lines  of  this  geometry 
is  shown  in  Fig.  4.  The  total  electrical  length  of  the 
line  is  1 14.6  mm  (4.51  in).  The  length  of  each  section 
is  16.6  mm  (0,65  in).  Figure  8  shows  the  signal  at  the 
receiver  end.  Comparison  is  made  to  the  control  line. 
(In  all  the  experiments  discussed  in  this  paper,  the 
control  line  is  a  straight  line  with  electrical  length 
equivalent  to  that  of  the  delay  line  under 
examination.  The  length  of  a  delay  line  is  considered 
to  be  the  length  of  a  line  that  runs  through  the  center 
of  the  line,  including  bends.) 

Observing  Fig.  8,  we  notice  that  the  wave 
tracing  analysis  predicted  the  ladders  that  precede  the 
main  signal  only.  While  such  prediction  is  quite 
important,  it  did^not  account  for  the  unexpected  high 
signal  level  received  after  the  arrival  of  the  main 
signal  (overshoot).  Although  the  high  level  occurs 
beyond  the  logic-switching  instant,  nevertheless,  it 
can  have  an  adverse  effect  when  -the  signal  switches 
to  zero  as  false  logic  could  be  triggered. 
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Fig.  7.  Wave  form  at  the  receiving  end  of  a 
serpentine  line  due  to  a  pulse  with  a  duration  of  T.  (a) 
T>2td.  (b)T<2td. 

For  the  second  test,  we  consider  a  more  realistic 
seipentine  line,  where  the  total  electrical  length  of  the 
line  is.  405.4  mm  (15.96  in).  Two  variations  of  this 
line  are  considered,  as  shown  in  Figs.  9  and  10,  and 
will  be  referred  to  as  Case  A  and  Case  B, 
respectively.  The  difference  between  the  two 
topologies  considered  is  that  one  has  shorter  sections 
than  the  other.  Notice  that  in  both  cases,  the  lines 
have  equivalent  length,  have  identical  line  separation 
between  adjacent  segments  and  both  topologies 
occupy  equal  board  area.  The  excitation  waveform  is 
a  pulse  having  a  finite  duration  of  1.4  nanoseconds 
and  rise  time  of  1 00  picoseconds. 

Figure  1 1  shows  the  response  of  these  two  lines 
in  comparison  to  the  reference  line.  The  serpentine 
line  designated  as  Case  B  has  longer  sections  and, 
consequently,  its  receiver  signal  had  ladders  that  are 


longer  than  Case  A.  However,  this  difference  is  of 
minor  importance  since  the  low-to-high  switching 
occurs  at  approximately  the  same  time. 


time  (nanosecond) 


Fig.  8.  Received  wave  form  for  the  seven-section 
serpentine  line. 


Fig.-9.  19-section  serpentine  line. 
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Fig.  10.  1 7-section  serpentine  line. 


0  1  2  3  4  9  B 

Ikne  (norrosecond) 


Fig.  1 1 .  Response  of  the  1 5.96  in  long  serpentine 
lines  in  comparison  to  the  reference  line. 

IIL  The  Spiritl  Delay  Line 

In  the  serpentine  line,  the  crosstalk  was  found  to 
accumulate  synchronously.  This  accumulation  can  be 
significant  enough  to  trigger  false  logic.  Ideally,  this 
crosstalk  needs  to  be  eliminated;  however,  since  the 
crosstalk  is  inversely  proportional  to  the  separation 
between  the  lines,  the  only  way  to  eliminate  or 
reduce  the  crosstalk  would  be  to  increase  the 
separation  between  the  lines  and/or  to  increase  the 
coupling  between  the  line  and  the  reference  plane. 
Increasing  the  separation  of  the  segments 
unfortunately  requires  larger  circuit  board  area, 
which  can  be  either  expensive,  or  impossible  in  light 
of  the  density  requirements.  Therefore,  given  that  for 


the  transmission  line  density  (or  etch  density)  to 
remain  unchanged,  the  separation  between  the  lines 
must  not  be  changed.  An  alternate  design,  which 
forces  the  crosstalk  to  accumulate  asynchronously, 
was  proposed  in  [3]  and  is  shown  in  Fig.  12.  This 
alternate  design,  the  flat  spiral  delay  line,  minimizes 
the  periodicity  in  the  transmission  line  routing  as 
much  as  practicable. 


A 


Fig.  1 2.  Spiral  delay  line. 

The  most  prominent  feature  of  the  spiral  delay 
line  is  the  spreading,  over  time,  of  the  crosstalk  noise. 
To  demonstrate  how  the  coupling  mechanism  works 
in  the  spiral  line,  we  will  initially  assume  that 
coupling  between  non-adjacent  lines  is  negligible. 
Let  us  consider  the  line  shown  in  Fig.  12,  and  let  us 
assume  that  the  signal  arrives  at  point  A  at  time  t=0. 
Once  the  signal  arrives  at  A,  forward  propagating 
crosstalk,  is  induced  at  point  B.  This  crosstalk 
arrives  before  the  main  signal  by  the  time  delay" 
needed  to  propagate  through  the  entire  length  of  the 
line  minus  the  segment  L26.  In  other  words,  the 
induced  forward  propagating  crosstalk  at  B  arrives 
only  after  the  time  needed  to  travel  the  segment  L26. 
(Henceforth,  the  time  delay  through  segment  LN  will 
be  denoted  by  ttw).  So  the  first  crosstalk  induced  at  B 
by  the  arrival  of  the  signal  at  A  arrives  at  t=tL26* 
When  the  main  signal  arrives  at  point  C,  forward 
propagating  crosstalk  is  induced  at  D,  arriving  at  the 
receiver  at  approximately  t==tLi+tL25+tL26* 

Similarly,  when  the  main  signal  arrives  at  E,  it 
induces  at  point  F  a  forward  crosstalk  which  arrives 
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at  the  receiver  at  t=  tL2+  tL26+  tL24-  Next,  the 

signal  arrives  at  G,  inducing  forward  propagating 
crosstalk  at  H,  which  arrives  in  t=  tL2+  tL3+  tL26+ 
tL25+  tu4  +  tL23.  Howcvcr,  the  signal  arriving  at  H 
will  also  induce  forward  propagating  crosstalk  at  L 
The  crosstalk  induced  at  I  arrives  at  the  receiver  at  t= 
tLi+  tL2+  tL3.  As  the  signal  propagates  further  towards 
the  receiver,  it  continues  to  induce  crosstalk  in  the 
adjacent  lines. 

The  above  wave-tracing  analysis  ignores 
coupling  between  non-adjacent  lines.  However,  even 
if  this  coupling  is  included,  it  is  easy  to  see  that  the 
crosstalk  noise  will  be  distributive  and  will  not 
accumulate  synchronously.  When  the  voltage 
excitation  is  a  pulse  of  finite  duration,  the  opposite 
polarity  crosstalk  induced  by  the  trailing  edge  of  the 
excitation  pulse  will  help  in  reducing  the  final 
cumulative  crosstalk. 

To  demonstrate  performance  of  the  flat  spiral 
delay  line,  we  construct  a  line  with  electrical  length, 
line  separation,  and  total  board  area,  all  equal  to  the 
serpentine  lines  shown  in  Figs.  9  and  10.  The  cross- 
section  parameters  are  given  in  Fig.  4,  and  the  rise 
time  and  pulse  duration  are  as  before  (pulse  width  of 
1.4  nanoseconds  and  rise  time  of  100  picoseconds). 
The  dimensions  of  this  spiral  line  are  shown  in  Fig.  - 
12. 

Figure  13  shows  the  signal  at  the  receiving  end 
of  the  spiral  line  in  comparison  with  the  control  line. 
The  outstanding  performance  of  the  spiral  line  is 
clearly  visible.  We  observe  that  the  crosstalk  noise 
was  spread  over  time  ahead  of  the  main  signal 
resulting  in  a  high  fidelity  signal.  In  fact,  for  the  case 
considered,  we  notice  ^at  the  maximum  crosstalk 
stays  at  or  below  10%  of  the  signal  level,  thus 
considerably  reducing  the  potential  for  triggering 
false  logic.  We  note  here  that  slight  shift  between  the 
signal  due  to  the  control  line  and  the  unadulterated 
signal  due  to  the  spiral  line  is  due  to  -a  slighr^ 
difference  in  line  lengths. 

The  singular  feature  of  the  spiral  line  is  its 
ability  to  spread.the  noise  over  a  larger  duration.  A 
key  observation  is  that  the  signal  arriving  at  the 
receiving  end  is  composed  of  the  unadulterated  pulse 
in  addition  to  the  asynchronous  noise,  This  implies 
that  the  spiral  line  can  lead  to  an  almost  perfect 
isolation  of  the  accumulative  noise  from  the 
unadulterated  pulse. 


Fig.  13.  Response  of  the  15.96  in  long  spiral  line  in 
comparison  to  the  reference  line. 


IV.  The  Concentric  Cs  Delay  Line 

From  the  analysis  and  simulation  of  the  spiral 
line,  the  concept  of  asynchronous  coupling  is 
emerging  as  a  powerful  mechanism  to  reduce  the 
apparent  dispersion  of  a  delay  line.  This  section 
introduces  a  novel  design  that  exploits  this 
mechanism  of  asynchronous  coupling.  The  new  delay 
line,  shown  in  Fig.  14,  will  be  referred  to  as  the 
concentric  Cs  delay  line.  (The  name  concentric  Cs 
was  chosen  since  this  new  transmission  line 
resembles  different  lines  shaped  as  the  letter  C  and 
centered  at  one  location.)  The  design  of  the  new  line 
was  motivated  by  the  theme  of  reducing  the 
periodicity  in  the  topology  of  the  structure  in  order  to 
minimize  the  accumulative  or  synchronous  coupling. 

Using  wave-tracing  analysis,  one  can 
qualitatively  predict  the  performance  of  the  Cs  line. 
Similar  to  the  spiral  line,  the  Cs  line  distributes  the 
coupling  over  time.  Consider  a  pulse  entering  the  line 
at  die  beginning  of  section  LI  (see  Fig.  J4).  This 
-^ulse  will  induce  a  crosstalk  that  will  prec^e  the 
primaiy  (or  unadulterated)  signal  by  a  time  period 
corresponding  to  the  propagation  over  segments  LI, 
L2,  L3,  L4,  L5  and  L6.  Notice  that  in  contrast  to  the 
spiral  line,  this  (first-appearing)  crosstalk  will 
precede  the  unadulterated  signal  by  a  relatively  short 
period,  however,  the  topology  of  the  Cs  line  prevents 
synchronous  accumulation.  One  can  qualitatively 
describe  the  periodicity  in  the  Cs  line  structure  as 
lying  in  between  that  of  the  serpentine  line  and  the 
spiral  line.  The  relative  advantage  of  the  mid-level 
periodicity  in  the  structure  becomes  more  apparent 
when  observing  the  full-wave  behavior  as  shown 
below. 
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To  demonstrate  the  behavior  of  the  Cs  line,  we 
design  a  line,  equivalent  in  electrical  length,  line 
spacing  and  total  board  area,  to  the  spiral  and 
serpentine  lines  studied  earlier.  The  dimensions  of 
this  design  are  shown  in  Fig.  14.  Figure  15  presents 
the  result  of  the  FDTD  simulation,  under  previous 
excitation  conditions,  showing  the  signal  at  the 
receiver  end. 

Ro 


L2 


1212  mm 


L5 


14 

L3 


Fig.  14.  The  concentric  Cs  delay  line. 

Comparison  is  made  to  the  spiral,  serpentine 
and  reference  lines.  Clearly  visible  from  Fig.  15  is  the 
concentration,  or  accumulation  of  the  noise  in  the 
close  proximity  of  the  unadulterated  signal.  We 
further  observe  that,  for  the  particular  topology  and 
line  dimensions  considered,  that  the  receiver  signal 
remains  at  or  below  35%  of  the  peak  amplitude  of  the 
transmitted  pulse. 


Fig.  15.  Response  of  the  concentric  Cs  line  in 
comparison  to  the  serpentine  and  spiral  lines  of  equal 
length. 


In  comparison  to  the  spiral  line,  the  Cs  line 
keeps  the  noise  closer  to  the  main  signal,  whereas  the 
spiral  line  distributes  the  noise  over  longer  time 
duration.  From  this  perspective,  the  spiral  line  is 
superior  to  the  Cs  line.  However,  if  one  were  to 
consider  an  excitation  composed  of  a  train  of  pulses, 
as  would  be  the  case  in  practical  scenarios,  it  might 
be  more  advantageous  to  force  the  crosstalk  noise  to 
accumulate  closer  to  the  pulse  that  generates  it  and 
not  interfere  with  other  pulses. 

V.  Summary 

This  paper  presented  a  qualitative  analysis  of 
serpentine  and  flat  spiral  delay  lines  based  on  the 
simply,  yet  powerful,  ray  tracing  technique.  The 
three-dimensional  FDTD  method  was  used  to  predict 
the  foil-wave  performance  of  these  delay  lines,  thus 
accounting  for  higher-order  modes,  multiple  and 
feedback  coupling  and  the  effect  of  orthogonal 
segments.  Despite  the  strength  and  completeness  of 
the  FDTD  method,  the  ray  tracing  method  facilitated 
an  understanding  of  the  coupling  mechanism  in 
complex-shaped  delay  lines  and  lead  to  the 
introduction  of  novel  delay  line  such  as  the  flat  spiral 
line.  A  new  line,  the  concentric  Cs  delay  line  was 
introduced  based  on  the  concept  of  reducing  the 
periodicity  in  the  structural  topology,  thus  forcing  the 
crosstalk  noise  to  accumulate  asynchronously. 
Numerical  simulation  using  the  three-dimensional 
foil-wave  FDTD  method  showed  that  the  new 
designs  result  in  receiver  waveform  that  is  less 
susceptible  to  triggering  false  logic  than  in  the  case  of 
the  serpentine  line. 
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An  Empirical  Approach  for  Design  of 
Wideband,  Probe-Fed,  U-Slot  Microstrip 
Patch  Antennas  on  Single-layer,  Infinite, 
Grounded  Substrates 

V.  Natarajan  and  D.  Chatterjee 


/lAj/rffcr— Wideband  mIcrostrIp  antennas  with  relatively 
simple  topologies  continue  to  attract  attention  for  design  of 
compact,  high-pcrformancG  communication  systems.  The 
coaxiallyTcd,  rectangular  patch  U-sIot  has  recently  been 
investigated  numerically  and  experimentally,  and  shown  to 
yield  10  dB  return>loss  bandwidths  in  excess  of  20%.  How¬ 
ever  there  are  no  analytical  models,  nor  any  systematic  de¬ 
sign  procedures  currently  available  that  can  aid  realizing 
these  conBgurations.  To  that  end,  based  on  extensive  CAD 
simulation  results  for  a  wide  range'^of  commercially  avail¬ 
able  microwave  substrates  (cr  =  2.94  to  10.2),  an  empirical 
design  methodology  is  derived  and  illustrated  by  examples. 
It  is  shown  that  the  present  empirical  design  technique, 
with  its  attendant  limitations,  generate  wideband  U-Slot 
designs  that  are  optimized  using  CAD  tools  such  as  IE3D 
within  a  few  iterations,  resulting  in  substantially  reduced 
overall  process  cycles. 

I.  INTRODUCTION 

The  theory  and  design  of  a  wide  variety  of  probe- 
fed,  microstrip  patch  antennas  for  various  applications 
has  been  well  documented  [1].  For  portable  phone  sys¬ 
tems  there  is  a  need  for  low-profile  (embedded)  antennas 
with  a  10  dB  return-loss  bandwidth  >  10%  in  addition 
to  other  desirable  electrical  characteristics  [2,  p.  3 1 2]. 
(The  return  loss,  for  this  paper,  is  taken  as  -201ogi(,  |r| 
in  dB,  where  P  is  the  reflection  coefficient.)  As  found 
in  [I,  ch.  9],  studies  on  wideband  microstrip  antenna 
designs  emphasize  techniques  such  as  multi-layer  sub¬ 
strates,  parasitic  elements  and  apcrture-ccwpled  excita¬ 
tions.  However  such  approaches  obviate  the  realization 
of  low-profile,  compact  antenna  topologies,  or  may  com¬ 
plicate  the  fabrication  process  due  to  the  need  for  sophis¬ 
ticated  feed  element  design(s)  [3]. 

Interestingly  the  U-slot  antenna,  first  reported  in  [4], 
was  a  new  form  in  ultrawideband  microstrip  antenna  de¬ 
sign  since  it  could  generate  «  40%  bandwidth  by  main¬ 
taining  very  simple  feed  and  patch  designs  on  single- 
layer  foam  substrates.  For  wideband  applications  it  thus 
appears  that  the  U-slot  design  is  a  pioneering  concept  as 
it  is  indeed  a  very  formidable  alternative  to  the  existing 
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wideband  patch  topologies  [1],[3].  The  subject  of  this 
paper  is  to  further  explore  some  advancements  in  design 
of  wideband,  probe-fed,  U-slot  patches  on  single-layer 
substrates. 

Most  results  for  this  novel  design  are  available  for  air 
(c,.  =:  1)  or  foam  (Cr  »  1)  substrates  [4]-(7].  One  finds 
appropriate  results  for  Cr  =  2.33  in  [8],  inclusive  of  fi¬ 
nite  ground-plane  and  substrate  truncation  effects.  Cur¬ 
rently  there  are  no  analytical  models,  nor  empirical  de¬ 
sign  relations  available  to  initiate  an  U-sIot  design  from 
some  nominal  specifications.  In  [6]  an  attempt  has  been 
made  to  use  a  [f]3x3  matrix  representation  for  the  U- 
slot,  but  the  analytical  expressions  for  the  diagonal  el¬ 
ements  Yn  ,22,33  are  unavailable.  (It  appears  from  [6] 
that  their  determination  was  done  using  iterative  exper¬ 
imental  techniques).  Explicit  formulas  for  the  two  res¬ 
onant  frequencies  of  the  U-slot,  with  validation  results 
for  foam  substrate  are  available  in  [9],[10]  that  are  good 
only  as  important  checks  in  a  simulation. 

Reduced  U-slot  patch  size  topologies  have  been  re¬ 
ported  in  [1I]-[14]  with  an  average  bandwidth  of  ~ 
30%.  Furthermore,  10  dB  return  loss  bandwidths  of  , 
44%  and  50%  have  been  reported  in  [15]  and  [16],  re¬ 
spectively.  Dual-band  designs  have  been  reported  in  [1 7] 
for  substrates  with  permittivity  Cr  ^  4.4,  and  in  [18] 
wideband  patch  designs  on  multi-layer  substrates  have 
been  reported.  Results  for  U-Slot  performance  on  mi¬ 
crowave  substrates  have  been  summarized  in  [19],  but 
without  any  design  information.  Recently,  a  different 
design  procedure  for  U-SIot  on  single-layer  microwave 
substrates  has  been  reported  in  [20].  However,  there  are 
significant  differences  between  [20],  and  the  proposed 
approach  in  this  paper.  These  differences  will  be  identi¬ 
fied  later  in  this  paper. 

The  information  gleaned  from  [21]  suggests  avoiding 
use  of  moment-method  [22], [23]  based  CAD  tools  like 
IE3D  [24]  for  initiating  single-element  patch  designs, 
due  to  prohibitively  high  computational  cost.  Thus  effi¬ 
cient  design  processes,  despite  their  heuristic/empirical 
nature,  can  still  help  the  overall  simulation  cycle  be  cost- 
effective.  Normally  such  empirical  procedures  reduce 
substantial  savings  in  computational  resources  by  requir¬ 
ing  fewer  iterations  in  the  final  CAD  optimization  of  the 
antenna  topology  [21].  Since  no  formal,  systematic  pro- 
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cedures  are  currently  available  for  the  design  of  U-SIot, 
the  purpose  of  this  investigation  is  to  develop  guidelines 
and  present  empirical  formulas  to  aid  in  realizing  such 
goals. 

The  empirical  formulas  developed  in  this  paper  apply 
to  probe'fed  designs  on  single-layer  grounded  substrates 
that  are  (ideally)  infinite  in  extent.  In  addition,  empiri¬ 
cal  design  formulas,  obtained  via  moment-method  based 
parametric  simulations,  apply  to  U-slots  patch  designs 
with  definite  geometrical  symmetry  as  elaborated  later  in 
this  paper.  The  contents  of  the  paper  are  outlined  next. 

To  that  end,  following  [25]-[27],  validation  results  for 
IE3D  against  appropriate  measured  data  for  microwave 
substrates  (Cr  =  2.33)  from  [8],  are  included.  The 
IE3D  code  validation  results  are  followed  by  a  care¬ 
ful  analysis  of  the  various  U-Slot  designs  studied  earlier 
[8],[28],  resulting  in  various  dimensional  invariance  rela¬ 
tionships  that  are  crucial  for  U-Slot  design  on  microwave 
(and  foam/air)  substrates.  Selected  results  demonstrat¬ 
ing  the  validity  of  the  empirical  formulas  are  included 
from  [29], [30].  Finally  the  major  observations  are  sum¬ 
marized,  and  a  list  of  relevant  references  is  included. 

II.  Problem  Description 

In  Fig.  1  the  geometry  of  a  probe-fed,  U-Slot  patch  on 
a  single-layer  substrate  is  shown  with  all  the  dimensions 
indicated  therein.  This  topology  is  a  simple  modification 
to  a  probe-fed  rectangular  patch  antenna,  the  latter  being 
generally  a  narrowband  radiating  element  [1]. 


Fig.  1.  Physical  topology  of  a  coaxially-fed,  single-layer,  rectangular 
patch  U-Slot  microstrip  antenna 

In  absence  of  any  analytical,  i.e.,  cavity  or  transmis¬ 
sion  line  models  (1,  chs.  4,5],  one  can  still  investi¬ 
gate  the  effects  on  a  performance  characteristic  (such 
as  impedance,  gain,  etc,)  due  to  variations  in  sub¬ 
strate/patch  geometry  via  careful  measurements  or  rig¬ 
orous,  full-wave  CAD  simulations  [21].  Development  of 
rapidly  iterative  design  procedures  could  involve  heuris¬ 
tic/empirical  approaches,  subject  to  further  refinements 
via  optimizations  [21]. 


For  the  investigation  reported  here,  the  main  aim  is 
to  examine  how  parameters  such  as  substrate  thickness, 
overall  patch  dimensions,  slot  width,  probe  location  and 
radius,  as  shown  in  Fig.  1,  affect  the  wideband  perfor¬ 
mance.  The  generic  nature  of  the  impedance  character- 


Fig.  2.  Typical  impedance  loci  characteristics.  The  perfonnance  of 
any  wideband  design  (or  modifications)  is  desired  such  that  the  loop 
of  the  impedance  loci  1,  2  &  3  encircle  the  center  (VSWR  =  1)  of  the 
Smith  Chart  as  in  locus  4. 

istics  is  shown  via  a  Smith  Chart  in  Fig.  2.  The  desired 
characteristic  is  depicted  in  #  4.  The  wideband  behavior 
of  the  antenna  will  be  superior  if  the  size  of  the  loop  for 
the  impedance  locus  #  4  shrinks  to  the  VSWR  =  1  point 
on  the  Smith  Chart.  (In  addition,  most  of  the  frequencies 
of  interest  has  to  lie  on  that  loop.)  For  practical  applica¬ 
tions,  the  size  and  location  of  loop  in  an  impedance 
loci  should  be  such  that  the  VSWR  <  2,  corre^onding  to 
a  return  loss  of  10  dB. 

Generally,  the  impedance  loci  will  be  far  removed 
from  the  desired  behavior  shown  in  #  4,  re.,  it  will  be 
more  like  the  loci  #  1,  2  or  3.  The  wideband  problem 
then  reduces  to  the  study  of  how  the  changes  in  various 
dimensions  in  Fig.  1  could  transform  loci  #  1,  2  &  3 
such  that  a  loop  could  be-obtained  in  the  impedance  loci 
meeting  the  criterion  VSWR  <  2. 

It  is  po'Ssible  to  analyze  the  impedance  behavior  using 
analytical  (cavity  model)  expressions,  to  a  very  good  de¬ 
gree  of  accuracy  [1,  chs.  4  and  5].  Such  would  facilitate 
rapid  parametric  simulations,  prior  to  any  computation¬ 
ally  intensive,  full-wave  analysis.  Since  there  exists  no 
such  analytical  model  for  U-Slots,  having  recourse  to  an 
alternate  route  for  rapid  parametric  studies  appears  crit¬ 
ical  before  any  CAD-based  optimization.  To  that  end,  it 
is  necessary  to  examine  the  capability  of  the  IE3D  code 
for  numerical  modeling  of  U-Slot  geometries.  The  re¬ 
sults  for  the  IE3D  code  validation  are  shown  next. 

III.  Validation  Results  for  the  ie3D  Code  [24] 

As  mentioned  in  [25]-[27],  a  CAD  tool  needs  to  be 
validated  for  an  ensemble  of  appropriate  test  cases  that 
are  closest  to  the  topology  being  studied.  Furthermore, 
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since  code  validation  is  an  open-ended  process,  a  judi¬ 
cious  selection  of  the  test  cases  forms  an  important  part 
of  such  investigation.  To  that  end,  it  was  decided  to  ex¬ 
amine  the  capabilities  of  the  IE3D  code  against  the  mea¬ 
sured  radiation  pattern  data  (in  <l>  =  0°^90°  planes)  for 
U-Slots  fabricated  on  substrate  €r  =  2.33  as  given  in  [8, 
Fig.  5]  corresponding  to  a  frequency  of  3.56  GHz.  To 
the  best  of  the  knowledge  of  the  present  investigators, 
reference  [8]  is  the  only  source  for  which  measured  and 
computed  data  are  available  for  U-Slots  on  finite  ground 
planes  for  microwave  substrates.  (Most  of  the  data,  as 
mentioned  earlier,  is  for  foam  (or  air)  substrates  for  U- 
Slot  topologies.) 

The  dimensions  of  the  antenna  B  as  in  [8,  Table  I] 
arc:  W  =  3.6cms,L  =  2.6cms,W8  =  l.dcmSjL#  = 
1.8cms,b  =  0-4cms,t  =  0.2cms,Xp  =  0.0cms,Yp  = 
0.0  ems,  and  F  =  1 .3  ems,  referring  to  Fig.  1  herein.  The 
radius  of  the  probe  could  not  be  found  in  [8];  after  sev¬ 
eral  trials,  it  was  found  that  dp^obe  =  0.127  ems  in  Fig. 
1  provided  the  best  agreement  with  the  data  in  [8].  The 
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Fig.  3.  Measured  -  ©  -  ©  -  data  [8.  Fig.  5];  computed  -  O- 
{IE3D  8.0)  and  -  *  -  •  -  (1 E3D  9.0),  for  total  fields  in  0  =  0®  plane 

IE3D  results  in  Figs.  3  and  4  arc  for  a  12  x  10  ems  rect¬ 
angular  ground  plane.  The  agreement  between  measured 
[8]  and  simulated  results  are  reasonably  acceptable  at  all 
angular  regions  except  near  6  90^,270°.  The  rea- 

son(s)  for  the  discrepancies  are  explained  below. 

The  actual  topology  analyzed  in  [8,  Fig.  1]  was  a 
U-Slot  located  on  truncated,  rectangular  substrate  on  a 
finite  rectangular  ground  plane.  In  the  IE3D  simula¬ 
tions,  the  radiating  patch  was  located  on  an  infinite  sub¬ 
strate  backed  by  a  finite  rectangular  ground  plane  of  the 
same  dimensions.  This  model  cannot  account  for  the 
surface  wave  diffraction  by  the  truncated  dielectric  sub¬ 
strate.  Since  surface  waves  arc  dominant  near  the  air- 
subslratc  interface,  (i.e.  9  90°,  270°),  the  radiation 

behavior  is  not  accurately  predicted  near  this  region  as 
seen  in  Figs.  3  and  4.  This  present  limitation  in  the 
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Fig.  4.  Measured  -  ©  -  ©  -  data  (8,  Fig.  5);  computed  -□  -  □- 
(IE3D  8.0)  and  ~  -  (IE3D  9.0),  total  fields  in  ^  =  90®  plane. 

IE3D  code  [24]  is  due  to  lack  of  implementation  of  the 
appropriate  microstrip  Green’s  function  that  can  account 
for  diffractions  due  to  dielectric  and  ground-plane  trun¬ 
cations.  In  [8],  since  the  FDTD  technique  was  used  to 
calculate  the  radiation  pattern,  such  edge  diffractions  are 
considered  in  situ.  Thus,  computed  patterns  in  [8]  do  not 
show  any  ficticious  discontinuities  near  9  90°,  270°, 

unlike  the  IE3D  results  presented  here. 

The  foregoing  results  demonstrate  the  limitations  of 
the  IE3D  code  when  applied  to  modeling  of  antennas 
on  grounded,  truncated  substrates.  However  the  2:1 
VSWR  bandwidth  of  antenna  B,  computed  via  IE3D, 
was  around  24%  -  in  good  agreement  with  measurements 
[8,  Table  II].  Also,  the  IE3D  code  had  been  used  to 
replicate  the  results  for  microstrip  antennas  on  infinite, 
grounded  substrates  with  the  test  cases  chosen  from  var¬ 
ious  topologies  in  [1].  In  all  these  cases  the  agreements 
were  very  good  with  published  data.  Since  the  scope  of 
this  present  investigation  is  limited  lo  infinite,  grounded 
substraies-the  type  of  discrepanciesin  Figs.  3  and  4  are 
not  likely  to  affect  the  results. 

IV.  DlMENSrONAL  INVARIANCE  IN  U-SLOT  DESIGN 

The  key  to  the  development  of  the  empirical  design 
procedure  is  the  establishment  of  the  dimensional  invari¬ 
ance  of  the  U-Slot  studied  in  [8]  and  [28].  These  results 
are  summarized  below  in  table  I  from  [29], [30].  In  table 
I  one  finds  that  the  only  parameter  which  changes  with 
substrate  is  and  all  other  dimensional  ratios  re¬ 
main  almost  invariant.  Consequently,  to  design  a  U-Slot 
on  an  infinite,  grounded  microwave  substrate  the  deter¬ 
mination  of  ^  for  a  specific  substrate  (Cr  and  h)  and 
resonant/design  frequency,  fr,  is  the  key  step.  One  can 
then  use  the  information  in  table  I  to  derive  the  topology 
of  the  patch  as  shown  in  Fig.  1 .  From  columns  3, 4  and 
7  in  table  I  one  can  easily  deduce  that  ^  w  1.38.  Inter- 
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TABLE  I 

Dimensional  invariance  in  U-Slot  designs. 


€r 

W"' 

T 

f: 

“51 — 

b 

- 2^ - 

w: 

— w — 

WT 

1,0 

8.168 

1.515 

0.835 

4.237 

0.13 

3.203 

2,33 

4.49 

1.445 

0.777 

4.5 

0.144 

2.573 

4.0 

3.87 

1.443 

0,776 

4.51 

0.144 

2.573 

9.8 

2.87 

1.442 

0.777 

4.48 

0.144 

2.574 

2.33 

5.624 

1.444 

0.777 

4.5 

0.143 

2.571 

The  data  have  been  obtained  from  | 

28]  and 

8],  for  the  first 

four  and  last  row,  respectively.  For  all  the  cases  cited  here, 
the  minimum  and  maximum  bandwidths  were  15%  and 
42%,  respectively. .  The  data  for  Cr  =  2.33  in  the  second 
and  fifth  rows  refer  to  U-Slot  topologies  from  [28]  and  [8], 
corresponding  to  900  MHz  and  3.26  GHz,  respectively. 


estingly,  this  fact  appears  to  have  been  confirmed  for  the 
U-Slot  data  presented  in  [1 9,  table  lb]. 

At  this  stage  it  is  important  to  distinguish  between  the 
approach  in  [20]  and  this  paper.  It  is  noted  that  [20],  like 
this  paper,  doesn’t  contain  any  full-wave  mathematical 
analysis  for  U-Slot.  One  of  the  main  differences,  in  con¬ 
text  of  table  I,  is  the  determination  of  U-SIot  dimensions 
[20,  sec.  III].  The  underlying  assumption  in  [20,  sec. 

Ill]  is  the  existence  of  four  different  resonant  frequen¬ 
cies  of  the  impedance  loop  (as  shown  in  locus  #  4  in  Fig. 

2)  for  the  U-Slot.  For  high  Cr  substrates  such  multiple 
resonances  may  not  occur,  but  an  impedance  loop  could 
still  form,  as  shown  in  Fig.2,  away  from  the  zero  reac¬ 
tance  (jX  =  0)  line  on  the  Smith  Chart.  Apparently,  this 
restricts  the  technique  in  [20,  sec.  HI]  primarily  to  low 
£r  substrates.  In  contrast,  the  dimensional  invariances 
(table  I)  apply  to  low,  medium  and  high  permittivity  sub¬ 
strates. 

Since  it  is  important  how  the  various  dimensions  in 
Fig.  1  could  affect  the  bandwidth,  a  detailed  study  was 
undertaken  to  examine  such  effects,  for  which  salient 
features  are  shown  in  the  following  section. 

V.  Parametric  Modeling  Studies  Via  ie3D  Code  [24] 

The  primary  objective  of  the  parametric  simulations 
is  to  examine  the  nature  of  the  input  impedance  varia¬ 
tions  as  shown  in  Fig.  2  in  section  II.  Assuming  that 
the  initial  topology  of  the  U-Slot  has  been  designed  us¬ 
ing  the  information  in  section  IV,  it  is  still  possible  that 
the  desired  bandwidth  may  not  have  been  achieved.  This 
implies  that  the  initial  design  needs  further  optimization, 
which  in  view  of  Fig.  2  implies  that  the  impedance  loop 
should  be  shrunk  to  encircle  the  vicinity  of  the  center  of 
the  Smid)  Chart.  The  parameters  that  exercise  significant 
control  on  the  impedance  loop  size  and  location  are  criti¬ 
cal  to  the  optimization  process.  Results  for  low  and  high 
permittivity  substrates  are  available  in  [29],  and  only  se¬ 
lected  results  for  Cr  =  4.5  are  shown  here  from  Figs. 

5  to  9.  The  data  are  included  in  the  individual  figures 
captions,  and  hence  are  not  repeated  here  to  avoid  te¬ 
dium.  In  Fig.  5,  Yp  =  -0.2  and  0.4  ems  refer  to  probe 
locations  below  and  above  the  origin  of  the  coordinate 


Fig.  5.  Effect  of  probe  location  on  the  impedance  behavior  of  U-Slol: 
€t  =  4.5,  tan  d  =  0.002,  W=  4.89,  L=  3.538,  h=  1.27,  L,  = 
2.45,  W,  =  1.9,  t=:  0.274,  a=  0.777,  b=  0.311,  =  2.5), 

<lpro6e  =  0.127  and  Xp  =  0.0  -  all  in  ems;  •  —  •  —  .  (Yp  =  0.4 
ems),  +  —  -b  -  +  (Yp  =  0.2  ems),  and  o  -  o  -  o  (Yp  =  -0.2  ems), 
referring  to  the  dimensions  shown  in  Fig.  1 . 


Fig.  6,  Effect  of  probe  radius  on  the  impedance  behavior  of  U-Slot: 
€r  =  4.5,  tan  <5  =  0.002,  W=  4,89.  L=  3.538,  has  1.27,  L,  = 
2.45,  W,  =  1.9,  t=  0.274,  a=  0.777.  b=  0,311,  (r.e.,f  =  2.5), 
Xp  =Yp  =  0.0  -  all  in  ems;  o  -  o  -  o  (dp^^ie  =  0.08636  emsX 

■  “  ‘ ’  (^ro6c  —  0.127  ems),  and  -| - b  —  +  (dprobe  —  0-2  ems), 

referring  to  the  dimensions  shown  in  Fig.  1 . 
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system,  respectively,  with  Xp  =  0,  as  shown  in  Fig.  1 . 
The  result  indicates  the  trend  that  as  the  probe  is  moved 
away  from  the  edge  of  the  slot,  the  impedance  loop  be¬ 
comes  more  inductive  and  its  size  decreases.  Similar 
trends  were  observed  for  other  substrate  cases  in  [29]. 

Fig.  6  shows  the  effects  of  the  probe  radius  (= 
|dprofir)  on  the  input  impedance.  In  contrast  to  the  re¬ 
sult  in  Fig.  5,  variations  in  probe  radius  doesn’t  shrink  or 
expand  the  size  of  the  loop.  The  comparison  further  indi¬ 
cates  that  the  dominant  effect  of  the  probe  on  the  U-Slot 
input  impedance  is  determined  by  its  location,  and  not 
radius.  Control  of  the  probe  radius  could  thus  be  viewed 
as  resulting  in  a  ‘fine-tuning’  mechanism  in  order  to  ob¬ 
tain  the  desired  wideband  behavior. 


1.0 


-1.0 


Fig.  7.  Effects  of  substrate  thickness  on  the  impedance  behavior  of 
U-Slot:  Cr  =  4.5,  tan  5  =  0.002.  W=  4.89.  L=  3.538,  L,  = 
2.45,  W,  =  1.9,  t=  0.274,  a=  0.777,  b=  0.311,  =  2.5), 

dprobe  =  0.127,  Xp  =Yp  =  0.0  .  all  in  ems;  o  -  o  -  o  (h=  0.6 

eras), - (h=  1.27  ems),  and  +  -  +  -  -Kh=  1.5  eras),  referring 

to  the  dimensions  shown  in  Fig.  1. 

In  Fig.  7  the  trends  in  impedance  behavior  with  in¬ 
crease  in  substrate  thickness,  h,  are  shown.  As  the  thick¬ 
ness  increases  from  0.6  ems  to  1 .5  ems,  the  impedance 
loop  decreases  in  size  and  becomes  more  capacitive  in 
character.  Forh  =  1.27  ems,  a  loop  in  the  impedance  be¬ 
havior  is  formed  closest  to  the  center  of  the  Smith  Chart 
-  indicative  of  wideband  behavior. 

Variation  in  the  U-slot  width,  t,  results  in  impedance 
changes  (Fig.  8)  similar  to  the  one  observed  for  probe- 
location  variations  in  Fig.  5.  As  the  U-slot  width  in¬ 
creases  from  0.2  ems  to  0.35  ems,  the  impedance  loop 
changes  from  being  inductive  to  capacitive,  and  for  a 
slot-width  t  =  0.274  ems  the  loop  is  located  close  to 
the  center  of  the  Smith  Chart. 

Increase  in  the  J  ratio  from  0.5  to  4.5  doesn’t  cause 
any  significant  change  in  the  location  of  the  impedance 
loop,  but  results  in  shrinking  of  its  size.  This  can  be 


1.0 


-1.0 


Fig.  8.  Effects  of  slot  width,  t,  on  the  impedance  behavior  of  U- 
Slot:  (r  =  4,5,  tan  J  =  0.002,  W=  4.89,  L=  3.538,  U  =  2.45. 
W,  =  1.9,  a=  0.777,  b=  0.311,  =  2.5),  dpro6e  =  0.127, 

Xp  =Yp  =5  0.0  -  all  in  ems; - (t=  0.2cms),  o-o-o(t=  0.274 

ems),  H - 1 - Kt=  0.35  eras)  referring  to  the  dimensions  shown  in 

Fig.  1. 


1.0 


-1.0 


Fig.  9.  Effects  of  r  ratio  on  the  impedance  behavior  of  U-Slot:  €r  = 
4.5,  tan  6  =  0.002,  W=:  4.89,  L=  3.538,L,  =  2.45,  Wj  =  1.9, 

t=  0,274,  dp^ie  =  0127,  Xp  =Yp  -  0.0  -  all  in  ems; - 

4.5),  o  -  o  -  o  (|  =  1.0),  H - 1-  “  +  (f  =  0.5),  referring  to 

the  dimensions  shown  in  Fig.  1. 
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inferred  from  Fig.  9. 

The  information  gleaned  from  the  most  important 
parametric  simulation  results,  as  shown  in  Figs.  5  to  9, 
suggests  the  following  optimization  guidelines  for  wide¬ 
band  U-Slot  design: 

(a)  change  the  slot  width,  t,  probe  location  Yp,  and 
substrate  thickness  h,  such  that  the  impedance 
loop  encircles  the  region  in  the  close  vicinity  of 
the  center  of  the  Smith  Chart; 

(b)  following  step  (a),  if  the  size  of  the  loop  is  unde¬ 
sirably  large  or  small,  increase  or  decrease  the  ^ 
ratio  to  reduce  the  loop  size  without  affecting  the 
location  of  the  loop  to  achieve  larger  bandwidth 

(c)  one  may,  optionally,  change  the  probe  radius  to 
move  the  impedance  loop  such  that  it  encircles,  or 
is  close  to  the  Smith  center  (VSWR  =  1) 

The  preceding  simulation  results  were  obtained  for  con¬ 
figurations  where  the  U-slot  is  symmetric  about  the  x 
axis  and  the  probe  is  located  such  that  Xp  =  0,  as 
shown  in  Fig.  1.  Again,  it  is  important  to  distinguish 
between  [20]  in  context  of  the  parametric  simulations  in 
the  present  investigation.  This  will  be  followed  by  the 
last  part,  i.e.,  development  of  the  empirical  design  equa¬ 
tions,  in  section  VI. 

The  parametric  simulations  in  [20,  sec.  11]  focussed 
mainly  on  the  variation  in  resonant  frequencies.  (The  U- 
Slot  geometry  studied  in  [4]  was  for  air  =  1  which 
was  scaled  to  Cr  =  2.2  in  [20]).  The  technique  in  [20, 
sec.  ni]  is  based  on  the  availability  of  limited  data.  Fur¬ 
thermore,  as  stated  in  [20,  sec.  the  effects  of  sub¬ 
strate  and  (probe)  feed  were  not  investigated  in  detail. 
The  information  gleaned  from  the  overall  comparisons 
between  [20]  and  the  present  paper  suggests  that  the  re¬ 
sults  included  here  have  broader  scope  of  applicability 
compared  to  [20]. 

VI.  Development  of  Empirical  Design  Formulas 

In  view  of  the  observations  on  dimensional  invariance, 
as  presented  in  table  I,  the  important  factor  that  initiates 
the  U-slot  design  is  a  knowledge  of  the  ^  ratio  fi'ora  an 
a-priori  knowledge  of  some  nominal  specifications. 

Consequently,  it  was  decided  to  examine  the  relation¬ 
ships  between  resonant  frequency,  fr,  substrate  parame¬ 
ters  €r  and  h,  and,  the  larger  dimension,  W,  of  theU-Slot. 
To  that  end,  following  the  data  in  table  I,  it  was  decided 
to  vary  ^  between  2  and  7,  in  increments  of  0.5.  A  typ¬ 
ical  substrate  was  chosen  and  from  the  pre-selected  ~ 
values,  the  U-Slot  dimensions  were  found  with  the  aid  of 
table  I.  Following  this  procedure,  various  U-Slots  were 
designed  for  a  wide  class  of  practical  substrates  avail¬ 
able  from  Rogers  Corp.  (We  must  emphasize  that  at  this 
stage  the  resonant  frequency,  /r,  is  unknown  and  was 
determined  as  described  below.) 

These  U-Slot  topologies  were  then  characterized  by 
the  full-wave  CAD  tool  IE3D  [24].  The  resonant  fre¬ 
quency  /r,  defined  by  zero  reactance  on  the  Smith  Chart, 
and  the  corresponding  fractional  2:1  VSWR  bandwidths 


for  each  case  were  noted,  (If  there  were  several  reso¬ 
nant  frequencies  in  the  2: 1  VSWR  range,  an  average  es¬ 
timate  of  fr  was  taken  [29].)  Each  discrete  pair  of  ^ 
and  fr  values,  corresponding  to  an  individual  design, 
were  plotted.  The  MATLAB  software  (version  6.1,  re¬ 
lease  1 2. 1 )  was  used  to  obtain  a  quadratic  relation  (‘best 
fit’)  for  these  ^  vs.  fr  plots.  For  a  any  specific  Cr, 
and  various  h,  several  such  equations  were  obtained,  as 
shown  in  table  II.  In  this  process,  it  was  observed  that 
a  >  20%  fractional  bandwidth  for  the  2:1  VSWR  range 
was  obtained  for  those  designs  obeying  3.5  <  ^  <  5.5. 
The  same  phenomenon  also  corresponded  to  the  range 
0.13  <  <  0.18  for  the  U-Slot  topologies  designed 

and  simulated  on  er  =  2.94, 4.5,  and  10.2,  (Here  A  cor¬ 
responded  to  the  resonant  frequency  ,  determined  from 
the  Smith  Chart.)  The  details  can  be  found  in  [29],[30] 
and  are  omitted  here  for  brevity. 

TABLE n 

Empirical  (quadratic)  equations  for  design  of  u-Slot 


■a 

0.635 

BiSi 

1.0 

f  =  0.2/* 
-2.4/r  +  9.2 

f  =  0.19/* 
-2.1/r  +  7.8 

f  =  -0.78// 
-0.3/.. +  8.8 

1.216 

f  =  0-62/? 

-4.8/r  +  12 

f  =  0-085// 
-1.7/r-f6.8 

f  =  -0.21// 
-2.2/.. +  8.1 

1.80 

-8.6/r  +  13 

^  =  0.89// 
-4.8/r-f8.5 

T^-b.89// 
+0.11/, +  4.5 

Here  fr  is  the  design  resonant  frequency  in  GHz,  and  h  and  e, 
are  the  substrate  thickness  and  permittivites,  respectively. 


The  next  section  illustrates  the  complete  empirical  de¬ 
sign  procedure  with  simulation  results  for  a  topology 
(Cr  =  3.27)  for  which  the  empirical  equations  are  not 
available  in  table  11. 

VII.  Empirical  Design  Technique  FoRjJ-SLOT 
In  this  sectionra  systematic  empirical  design  proce¬ 
dure  for  design  of  U-Slot  patch  antennas  on  microwave 
substrates  is  presented  from  [29].  It  is  shown,  from  the 
results  in  secs.  IV,  V  and  VI,  that  U-SIot  patch  an¬ 
tennas  can  be  realized  which  are  further  optimized  us¬ 
ing  IE3D  CAD  software  [24].  VSWR  and  boresi^t 
(^  =  O°,0  =  0®)  Gain  results  for  the  unoptimized  znd 
optimized\5-S\oX  topologies  are  included  to  demonstrate 
the  efficacy  of  the  empirical  design  procedure. 

The  limitation  of  this  design  procedure,  as  mentioned 
before,  is  that  the  U-Slot,  as  shown  in  Fig.  1,  is  located 
symmetrically  w.r.t  coordinate  axes  with  the  probe  is  on 
the  y  axis,  ^or  rapid  automated  calculations,  the  de¬ 
sign  procedure  can  easily  be  adapted  within  a  computer 
design/simulation  code.)  It  is  assumed  that  the  U-slot 
antenna  should  have  a  10  dB  return  loss  bandwidth  of 
>  20%,  after  final  optimization.  Another  limitation  is 


196 


197 


ACES  JOURNAL,  VOL.  1 8,  NO.  3,  NOVEMBER  2003 


that  the  dimension  a  =  6  in  this  design  approach,  and 
the  slot  width,  t,  remains  uniform. 

(1)  From  the  nominal  a-priori  specifications  for  reso¬ 

nant  frequency  fr,  select  a  commercially  available 
substrate  with  thickness  h  to  satisfy  the  crite¬ 
rion  0.13  <  <  0.18.  In  the  experience  of  the 

present  authors  the  upper  and  lower  limiting  values 
should  be  used  for  low-  and  high-permittivity  sub¬ 
strates,  respectively.  For  intermediate  permittivity 
substrates  (Cr  ^  ^-5)  the  criterion  0.15 

can  be  used. 

(2)  Employ  the  empirical  equations  from  table  II  to 
calculate  the  ^  ratio,  and  from  step  (1),  one  can 
subsequently  determine  the  overall  width  W.  (One 
may  check  for  the  additional  criterion  3.5  <  ^  ^ 
5.5,  upon  calculation  of  ^  ratio.) 

(3)  From  table  I,  one  uses  ^  «  2.57  to  determine 

Ws. 

(4)  From  table  I,  since  ^  w  0.777,  calculate  L,  with 
the  knowledge  of  T^s  from  step  (3). 

(5)  From  the  relation  ^  »  0.144  in  table  I,  calculate 
the  slot  width  t,  with  the  knowledge  of  W,  from 
(3).  (This  assumes  a  slot  of  uniform  width.) 

(6)  Similarly,  from  table  I,  via  the  relation  ^  «  4.5, 
calculate  b  with  a  knowledge  of  Ls  from  (4). 

(7)  Assume  a  =:  6  in  Fig.  1,  and  calculate  L  =  is  + 
(I  b  ^  Lig  4“  2q  “  is  “I"  26, 

(8)  Locate  the  coaxial  probe  exactly  at  the  center,  i.e., 
Xp  =  0  and  Yp  =  0,  (or  equivalently  F  =  f ). 

(9)  Simulate  the  U-Slot  geometry,  as  obtained  via 
steps  1  to  8,  using  IE3D  (or  any  other  microstrip 
CAD  package  [27]),  and  check  for  the  2:1  VSWR 
performance. 

(10)  By  examining  the  nature  of  the  impedance  vari¬ 
ation  on  a  Smith  Chart,  from  step  (9),  adjust  the 
parameters  (as  mentioned  in  sec,  V)  for  further 
improvement  of  wideband  performance  of  U-Slot. 

In  order  to  illustrate  the  application  of  the  preced¬ 
ing  design  steps,  several  topologies  were  modeled  -  for 
which  the  Jesuits  arc  available  in  [2^,[30].  The  results 
contained  in  [29]  (and  [30])'mainly  demonstrate  the  ap¬ 
plicability  of  the  empirical  design  procedure  for  those 
substrates  for  which  the  empirical  relations  are  contained 
in  table  II.  In  this  paper  separate  results  are  presented  for 
a  typical  case  fr  =  3.27  (TMM3),  for  which  no  informa¬ 
tion  is  available  in  tables  I  and  11.  Since  the  permittivity 
of  (TMM3)  lies  in  the  range  2.94  <  fr  =  3.27  <  4.5, 
so  either  of  the  equations  in  table  II  can  be  used.  These 
would  result  in  two  different  dimensions  for  the  U-Slot, 
as  ^  would  be  different  for  the  two  cases.  The  results 
from  IE3D  [24]  simulation  for  the  two  U-SIot  topologies 
arc  presented  here  to  illustrate  any  such  differences. 

For  the  TMM3  (Cr  =  3.27)  substrate,  an  operating  fre¬ 
quency  fr  =  2.3  GHz  was  chosen.  From  the  condition 
«  0.15,  it  was  found  that  h  =  1.08  cms.  Re¬ 
ferring  to  table  II  the  equations  for  fr  =  2.94  and  4.5, 
corresponding  to  a  substrate  thickness  h  =  1.0  cms. 


were  chosen.  It  was  found  that  ^  =  4.738(er  = 
2.94)  and  3.975 (fr  =  4.5)  via  the  two  appropriate  em¬ 
pirical  relations  from  table  11.  The  remainder  of  the  di¬ 
mensions  were  easily  found  following  steps  (1)  to  (8). 

At  this  stage,  the  (initial)  unoptimized  design,  ob¬ 
tained  by  following  steps  (1 )  to  (8)  only,  was  character¬ 
ized  by  the  IE3D  code.  The  VSWR  vs.  frequency  results 
were  then  examined,  and  the  probe  location  was  changed 
from  its  initial/unoptimized  value  (Xp  =  0  and  Yp  =  0) 
by  moving  the  probe  along  the  y  axis  to  Yp  =  0,1  cms, 
in  view  of  the  results  in  Fig,  5.  (Various  other  opti¬ 
mization  options  in  section  V  could  have  been  pursued 
as  well.)The  VSWR  and  boresight  gain  vs.  frequency 
for  the  unoptimized  and  optimized  U-Slot  topologies  are 
compared  to  demonstrate  the  quality  of  the  initial  (unop¬ 
timized)  design  obtained  via  steps  (1)  to  (8). 

The  final  dimensions  of  the  two  U-Slot  patches  (Cr  = 
3.27)  obtained  from  the  two  empirical  equations  in  ta¬ 
ble  II  are  given  below.  (These  geometries  include  the 
optimized  probe  locations  obtained  for  enhanced  band- 
widths.) 

(i)  via  the  empirical  equation  for  £,  =  2.94  in  table  II: 
W=  4.738,  L=  3.424,  W,  =  1.841,  L,  =  2.37, 
h=  1.0,  t=  0.2C5,  a=b=:  0.527,  Xp  =  0.  Yp  = 
Oil  anddpro6c  =  0.127  cms;  substrate Cr  =  3.27 

(ii)  via  the  empirical  equation  for  Cr  =  4.5  in  table  II: 
W=  3.975,  L=  2.872,  W,  =  1.545,  L,  =  1.988, 
h=  1.0,  t=  0.223,  a=b=  0.442,  Xp  =  0,  Yp  = 
0.1  and  dprobc  =  0.127  cms;  substrate  =  3.27 

For  the  two  topologies,  the  VSWR,  Gain  variations 
vs.  frequency  and  the  radiation  patterns  in  the  cardinal 
planes  (<l>  =  0°  and  90°)  were  obtained  via  IE3D  code 
[24],  and  are  shown  in  Figs.  10  to  1 7.  The  data  in  Figs. 
10  to  13  compare  the  performances  of  the  unoptimized 
and  optimized  designs.  The  results  are  briefly  discussed, 
next. 


Fig.  10.  lllustratirg  difTercnccs  in  VSWR  for  unoptimized  (*  —  •  —  »: 
Xp  =Yp  =  0,0  cms),  and  optimized  (o  —  o  —  o:  Xp  =  0.0  and 
Yp  =  0.1  cms)  U-SIot  geometry  as  described  in  (i). 


The  results  in  Figs.  10  and  1 1  suggest  that  the  un¬ 
optimized  U-Slot  design  for  Cr  =  3.27,  obtained  via 
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Fig.  1 1.  Illustiating  differences  in  boresight  Gain  (<^  =:  =  0^) 

for  unoptimized  (*  —  »  —  •;  Xp  =sYp  =  0.0  eras),  and  optimized 
(o  -  o  -  o:  Xp  =  0.0  and  Yp  =  0,1  eras)  U<S)ot  geometry  as 
described  in  (i). 

the  €r  =  2.94  equation,  performs  well  within  expecta¬ 
tions.  The  2:1  VSWR  bandwidths  for  the  unoptimized^ 
and  optimized  topologies  are  w  33%  and  30%,  respec¬ 
tively.  One  also  notices  that  for  the  unoptimized  case  the 
overall  VSWR  performance  is  somewhat  inferior  com¬ 
pared  to  the  opdmzied  U-Slot  topology.  This  conclu¬ 
sion  can  be  reached  by  examining  the  VSWR  behavior 
in  the  vicinity  of  2  GHz  for  the  two  cases.  By  moving  the 
probe,  however,  the  overall  VSWR  behavior  is  improved 
at  the  expense  of  some  reduction  in  bandwidth.  The  gain 
behavior  in  Fig.  1 1  does  not  show  marked  changes. 


Fig.  12.  Illustrating  differences  in  VSWR  for  unoptimized 
Xp  =Yp  =  0.0  ems),  and  optimized  (o  —  o  —  o:  Xp  =  0.0  and 
Yp  =  0.1  ems)  U-SIot  geometiy  as  described  in  (ii). 


Results  in  Figs.  12  and  13  exhibit  similar  characteris¬ 
tics  when  compared  to  Figs.  10  and  11.  Comparing  the 
VSWR  variations  in  Figs.  1 0  and  1 2,  one  notices  that  un¬ 
optimized  (initial)  design  is  a  good  estimate  that  could  be 
optimized  using  few  iterations  on  the  probe  location.  In¬ 
terestingly,  the  boresight  gain  variations  in  Figs.  1 1  and 
13  are  rattier  insensitive  to  probe  locations.  This  obser¬ 
vation  suggests  that  initial  (unoptimized)  U-SIot  designs. 


Fig.  13.  Illustrating  differences  in  boresight  Gain  =  0®,  =  0®) 
for  unoptimized  (*  —  *  —  *:  Xp  =Yp  =  0.0  ems),  and  optimized 
(o  —  0  —  o:  Xp  =  0.0  and  Yp  =  0.1  ems)  U-Slot  geometiy  as 
described  in  (ii). 

would  not  show  perceptible  changes  in  gain  behavior 
if  their  wideband  performance  is  enhanced  by  chang¬ 
ing  the  probe  location.  Similar  influence  of  probe  lo¬ 
cation  was  observed  for  other  designs,  and  are  contained 
in  [29], [30]. 


Fig.  14.  VSWR  characteristics  for  U-Slot  m  substrate  tr  = 
3.27, tan  6  =  0.002  and  h  —  1.0  cm;  ♦  —  *  —  *  and  o  —  o  —  o  refer 
to  desi^s  derived  from  empirical  equations  for  €r  =  2.94  and  4.5, 
respectively,  from  Table  II.  The  data  refers  to  optimized  designs  only. 

Results  in  Figs.  14  to  17  refer  to  qpriw/zei  topologies 
as  defined  earlier  in  (i)  and  (ii).  The  VSWR  results  in 
Fig.  14  indicate  that  both  designs  offer  «  30%  band- 
widths  corresponding  to  a  return  loss  of  1 0  dB.  However 
the  frequency  ranges  over  which  such  wideband  behav¬ 
ior  occurs  is  different  for  the  two  designs.  For  both  de¬ 
signs,  VSWR  <  2  at  the  original  design  frequency  of  2.3 
GHz. 

The  gain  characteristics  shown  in  Fig.  15  reveal  that 
the  two  U-Slot  designs  exhibit  overall  similarities,  except 
that  the  topology  derived  from  the  empirical  equation  for 
Cr  =  4.5  shows  improved ‘gain-flatness’.  However  at  2.3 
GHz  the  gain  values  are  0.6  and  4  dB,  respectively,  for 
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Fig.  15.  Boresight  (9  =  0®,^  =  0®)  gain  characteristics  forU-Slot 
on  sub.stratc  €r  =  3.27,  tan  <5  =  0.002  and  =  1.0  cm;  *  —  *  —  ♦ 
and  0—0  —  0  refer  to  designs  derived  from  empirical  equations  for 
(r  =  2.94  and  4.5,  respectively,  from  Table  11.  The  data  refer.?  to 
optimized  de.signs  only. 

U-SIo!  geometry  defined  in  (i)  and  (ii). 

Examining  the  VSWR  and  gain  characteristics  for  the 
two  optimized  designs  (i)  (tr  =  2.94),  and  (ii)  (Cr  = 
4.5),  one  can  explain  the  differences  by  realizing  that 
these  are  essentially  two  different  U-Slot  geometries  on 
the  same  substrate  material  (er  =  3.27  and  h  =  1.0  cms). 
When  these  two  were  simulated  via  the  IE3D  code  the 
results  were,  quite  predictably,  different  and  is  seen  in 
Figs.  Hand  15. 


0 


Fig.  16.  Radiation  pattern  vs.  polar  angle  9  in  ^  =  0“  plane  for 
U*Slot  on  substrate  (r  =  3.27,  land  =  0.002  and  /t  =  1.0  cm. 
0-0-0  (2.0  GHz)  and  □  -  □  -  □  (2.3  GHz)  refer  to  desi^s  derived 
from  empirical  equations  for  €t  =  2.94  and  4.5,  respectively,  from 
Table  II.  The  data  refers  to  optimized  designs  only. 

The  principal  plane  radiation  patterns  in  Figs.  16  and 
17.  The  patterns  arc  shown  for  2.0  and  2.3  GHz,  for  re¬ 
spective  maximum  boresight  gains  as  in  Fig.  15.  The 
results  show  almost  no  difference  for  the  two  U-Slot  de¬ 


signs. 
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Fig.  17.  Radiation  pattern  vs.  polar  angle  0  in  ^  =  90°  plane  for 
U-Slot  on  substrate  €r  —  3.27,  land  =  0.002  and  h  =  1.0  cm; 
0-0-0  (2.0  GHz)  and  □  —  a  —  □  (2.3  GHz)  refer  to  designs  derived 
from  empirical  equations  for  6r  =  2.94  and  4.5,  respectively,  from 
Table  n.  The  data  refers  to  optimized  designs  only. 

One  may  thus  conclude,  judging  qualitatively  the  re¬ 
sults  in  Figs.  10  to  17,  that  the  empirical  formulas  in 
table  IT  are  reasonably  reliable,  and  have  a  wider  scope 
of  applicability  as  compared  to  [20].  Follow-up  investi¬ 
gations,  comparing  the  present  design  method  and  [20, 
sec.  Ill]  are  planned  for  future  work. 

VIII.  Discussion  and  Future  Work 

The  empirical  design  technique  developed  here  for  U- 
slot  patches  in  section  VII  is  by  no  means  exhaustive. 
The  most  important  part  in  the  design  are  embodied  in 
steps  1  to  8  in  sec.  VII,  that  yield  a  basic  design.  Steps 
9  and  10  are  essentially  equivalent  to  the  optimization 
functions  available  in  the  IE3D  code  [24].  The  utility 
of  stepT9  and  10,  from  a  practical  point  of  view,  stem 
from  the  fact  that  not  all  microstrip  antenna  CAD  tools 
contain  this  optimization  facility  [27]  like  the  IE3D  code 
[24].  In  those  special  situations  steps  9  and  10'  can  play 
a  critical  role  in  the  final  design.  It  is  however  impor¬ 
tant  to  assess  how  the  optimizers,  currently  available  in 
IE3D,  compare  with  an  U-slot  design  obtained  via  the 
parametric  simulation  studies  (steps  9  &  1 0  in  sec.  VII). 
This  work  is  currently  in  progress,  and  the  final  (opti¬ 
mized)  results  shall  also  be  compared  against  the  An- 
soft  ENSEMBLE  microstrip  CAD  software  for  the  vari¬ 
ous  optimization  options  available  in  IE  3D  [24].  In  ad¬ 
dition,  in  view  of  the  recent  work  reported  in  [20],  ad¬ 
ditional  investigations  are  necessary  to  compare  the  two 
different  U-Slot  design  approaches,  i.e.,  in  section  VII 
and  [20,  sec.  HI].  This  comparison  needs  to  be  carried 
out  with  special  emphasis  on  high  permittivity  substrates 

fr  >  6.0. 
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As  mentioned  earlier,  this  empirical  design  technique 
is  limited  because  of  symmetries  associated  with  U-slot 
topologies.  While  the  dimensional  invariance  and  em¬ 
pirical  formulas  initiate  a  design,  the  complete  process 
doesn’t  provide  any  analytical  insight  into  the  electro¬ 
magnetic  behavior  of  the  antenna.  Recent  investigations 
[19], [20]  on  U-Slot  design  and  performance  modeling, 
suggest  the  need  for  analytical  development.  To  elimi¬ 
nate  this  present  limination,  efforts  to  develop  an  analyt¬ 
ical  formulation  -  based  on  the  generalized  cavity  [1 ,  pp. 
97-102]  and  multiport  network  [1,  pp.  103-108]  models 
-  are  under  consideration.  The  main  objective  of  such 
a  future  effort  would  be  to  develop/derive  semi-rigorous 
formulas  that  are  far  less  empirical  in  nature  than  that 
presented  here.  Such  effort  would  also  be  aimed  at  pro¬ 
viding  information  on  the  nature  of  currents  flowing  on 
the  patch  surface,  and  slot  resonanoe(s).  There  is  some 
information  on  slot  resonance  in  [9]  and  [10],  but  they 
appear  valid  for  low-permittivity  substrates. 

IX.  Conclusion 

Earlier  investigations  have  shown  that  introducing  an 
U-shaped  slot  on  the  radiating  surfece  of  a  probe-fed, 
rectangular  microstrip  patch  antenna,  on  a  single  layer 
substrate,  resulted  in  ultra-wideband  topologies  with  bet¬ 
ter  than  10  dB  return  loss  performance.  Since  most  of 
the  published  results  on  U-slot  were  for  foam  or  air  sub¬ 
strates,  and  no  systematic  design  procedure  was  avail¬ 
able,  a  technique  has  been  proposed  in  this  paper  for  the 
design  of  such  a  class  of  wideband  microstrip  antennas 
on  low  to  high  permittivity  microwave  substrates.  By 
examining  the  earlier  data  for  U-slots,  the  key  feature 
of  dimensional  invariance  was  established,  and  empiri¬ 
cal  equations  have  been  derived  based  on  data  obtained 
from  extensive  (IE3D)  simulations.  Subsequently,  a  sys¬ 
tematic  design  procedure,  for  rapid  parametric  simula¬ 
tion  and  design  of  U-slot  patch  antennas,  has  been  pro¬ 
posed  incoiporating  the  two  above-mentioned  features. 
Using  this  proposed  technique,  VSWR  and  boresight 

—  0°  and  9  =  0®)_gain  comparisons  were  performed 
for  the  unoptimized  and  optimized  U-Slot  topoFogiess-  - 
The  results  suggest  that  initial  designs  for  U-slots  that 
can  be  optimized  within  a  few  iterations  via  paramet¬ 
ric  simulations  or  state-of-art  microstrip  antenna  analysis 
CAD  tools,  such  as  IE3D  or  ENSEMBLE. 
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Abstract 

It  is  well  known  that  the  computation  of  magnetic  fields 
in  nonlinear  magnetic  media  may  be  carried  out  using 
various  techniques.  In  the  case  of  problems  involving 
complex  geometries  or  magnetic  media,  numerical 
approaches  become  especially  more  appealing.  The 
purpose  of  this  paper  is  to  present  an  automated  particle 
swarm  optimization  approach  using  which  field 
computations  may  be  carried  out,  via  energy 
minimization,  in  devices  involving  nonlinear  magnetic 
media.  The  approach  has  been  implemented  and 
simulations  were  carried  out  for  different  device 
configurations.  It  is  found  that  the  computations 
obtained  using  the  proposed  approach  are  in  good 
qualitative  and  quantitative  agreement  with  those 
obtained  using  the  finite-element  approach.  Details  of 
the  proposed  approach,  simulations  and  comparisons 
with  finite  element  results  are  given  in  the  paper. 


1.  INTRODUCTION 

It  is  known  that  magnetic  field  computation  in 
nonlinear  magnetic  media  may  be  carried  out  using 
various  techniques  (refer,  for  instance,  to  [1,2]). 
Obviously,  for  cases  involving  complex  geometries 
and/or  magnetic  media,  numerical  and  artificial 
intelligence  approaches  become  especially  more 
appealing  (see,  for  example,  [3,4]). 

Irrespective  of  the  adopted  approach,  geometrical 
domain  subdivision  is  usually  performed  and  local 
magnetic  quantities  are  considered.  One  way  to  obtain 
an  electromagnetic  field  solution  is  tluough  the 
minimization  of  the  problem’s  energy  functional,  which 
may  take  complicated  non-quadratic  forms.  The 
purpose  of  this  paper  is  to  present  an  automated  particle 
swarm  optimization  (PSO)  approach  [5,6]  using  which 
2-D  field  computations  may  be  carried  out  in  devices 
involving  nonlinear  magnetic  media.  More  specifically, 
the  magnetic  energy  functional  is  first  formulated  in 
terms  of  the  unknown  magnetic  vector  potentials 
corresponding  to  the  discretization  scheme. 

A  swarm  of  particles,  each  designated  by  a  position 
vector  that  represents  the  unknown  potentials,  is 
initially  randomly  generated.  These  position  vectors 
may  be  regarded  as  potential  solutions  to  the  energy 
minimization  problem.  The  swarm  is  repeatedly  moved 


(i.e.,  modified)  by  the  optimization  algorithm  and,  upon 
convergence,  the  unknown  magnetic  vector  potentials 
are  found.  Consequently,  the  field  distribution  is 
computed  everywhere. 

Among  the  advantages  of  the  proposed  approach  are 
its  simplicity,  ability  to  handle  complex  magnetic  media 
and  computational  efficiency.  The  proposed  approach 
has  been  implemented  and  computations  were  carried 
out  for  different  electromagnetic  device  configurations. 
These  computations  showed  good  qualitative  and 
quantitative  agreement  with  results  obtained  using  the 
finite-element  (FE)  approach.  Details  of  the  approach, 
computations  and  comparisons  with  results  obtained 
using  the  FE  approach  are  given  in  the  following 
sections. 


n.  PROBLEM  FORMULATION 

In  general,  2-D  electromagnetic  field  problems  may 
be  reduced  into  1-D  problems  using  the  following  well- 
known  magnetic  vector  potential  Au^  formulation: 

B=fJl={'7xAu^),  (1) 

where  fi  is  the  permeability,  jj^is  a  unit  vector 

orthogonal  to  the  problem  plane,  H  is  the  magnetic  field 
vector  and  B  is  the  magnetic  flux  density  vector. 

For  nonlinear  media,  the  magnetostatic  energy 
functional  E  may  be  expressed  in  the  form: 


!l^ 

-JA 
0  ^ 


ya 


(2) 


where  O  represents  the  problem  domain  area,  y  is  the 

reciprocal  of  the  magnetic  permeability  (i.e.,  Y= ), 
and  J  is  the  current  density  along  . 

Neglecting  hysteresis  effects,  the  B-H  relation  of 
most  non-linear  magnetic  materials,  especially  those 
used  in  electromagnetic  power  devices,  may  be 
reasonably  approximated  by: 


B  *  C  'J  if  eii  and  H  « 


C 

V  j 


(3) 
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where,  n  is  an  odd  number,  C  is  a  constant  while  eij 
and  ejjurc  unit  vectors  along  the  field  and  flux  density 
directions,  respectively. 

By  subdividing  the  problem  domain  into  finite 
magnetic  and  nonmagnetic  regions,  and  from 
expressions  (l)-(3),  the  magnetostatic  energy 
formulation  becomes: 


It  turns  out  that  obtaining  a  solution  of  a 
magnetostatic  field  problem  involving  non-linear 
magnetic  media  should  always  correspond  to  a 
minimization  of  expression  (7).  In  this  work,  obtaining 
the  solution  is  achieved  through  the  utilization  of  the 
PSO  approach  to  directly  search  for  the  appropriate 
vector  potential  values  as  explained  in  the  following 
section. 


r=/ 


f 

D 

db-JrAf'’ 

J 

0 

rft 

^  ^  J 

AQ;- 


in.  ENERGY  MINIMIZATION  USING 
PARTICLE  SWARM  OPTIMIZATION 


2fio 


-JrA^ 


AQr 


(4) 


where,  is  the  average  value  of  the  magnetic  vector 
potential  within  the  subdivision,  A^f.  is  the  area  of 
the  /'subdivision,  is  the  permeability  of  free  space, 
while  Pf„  and  represent  the  number  of  magnetic 
and  non-magnctic  domain  subdivisions. 

In  the  case  when  triangular  domain  subdivisions  are 
adopted,  the  value  of  the  vector  potential  within  the  r'* 
subdivision  may  obviously  be  expressed  by: 


Ar  (.X,  y)  =  OTto  (Arj.Ar2.Ar3  )+  (Arl.Ar2.Ar3  )* 

+ary(Ar].Ar2,Ar3)y,  (5) 

where  and  are  constants  within  the 

triangular  subdivision  and  may  be  formulated  in  terms 
of  the  vector  potential  values  at  its  three  vertices  (i.e., 


ArJ,Ar2>Ar3y 

Hence,  from  (1)  and  (5),  corresponding  flux  density 
components  can  be  written  in  the  form: 


^rx  +<^ryXArJ>Ar2.Af.3) 

Bfy  =  ~‘CXrx  (Afr] ,  Aj.2  ,  Af3  ) 

|^r|  “  » Ar2 » )+  (Af] ,  *  Aj-^ ) 


(6) 


Substituting  (6)  into  (4),  we  obtain: 


I' 

r=/ 


(Ar}>Af2^A^3  ^ry  (^r  7»^r2  AfJ  ^ _ 


-z 

r=] 


(w  +  2)C” 

-  /  Ur7  •^^r2 
'■  3 

^rxiArJ>Ar2fAfr3  )+  a^y  (A^},Af.2,Ar3 


r  +  Ar2  + 
-Jr  5 


Ail,.  (7) 


Particle  Swarm  Optimization  (PSO)  is  an 
evolutionary  computation  technique,  which  simulates 
the  social  behavior  of  insect  swarms  or  bird  flocks. 
Different  from  traditional  search  algorithms,  such  an 
evolutionary  computation  technique  works  on  a 
population  of  potential  solutions  of  the  search  space. 
Through  cooperation  and  competition  among  the 
potential  solutions,  this  technique  is  suitable  for  finding 
global  optima  in  complex  optimization  problems  [5]. 

The  idea  of  our  PSO  implementation  is  that  the 
whole  swarm  proceeds  in  the  direction  of  the  swarm 
member  with  the  best  fitness  in  a  more  or  less 
stochastic  way.  A  swarm  consists  of  several  particles 
Af,  where  each  particle  keeps  track  of  its  own  attributes. 
The  most  important  attribute  of  each  particle  is  its 
current  position  as  given  be  an  TV^-dimensional  vector  A* 

=  (Ay  ,A2v*-mA^)»  corresponding  to  a  potential 

solution  of  the  function  to  be  minimized. 

Along  with  the  position,  each  particle  has  a  current 

velocity,  v*=  (i/^.V2  >•— »  keeping  track  of  the 

speed  and  direction  in  which  the  particle  is  currently 
traveling.  Each  particle  also  has  a  current  fitness  value, 
which  is  obtained  by  evaluating  the  objective  function 
at  the  particle’s  current  position.  Additionally,  each 
particle  remembers  its  own  personal  best  position,  p* 

=  (Pj  y  P2 '  Pn^  ‘  swarm  level,  the  best 

overall  position  among  all  particles,  ,  is  also 
recorded.  The  index  of  the  best  particle  is  represented 
by  the  symbol  g.  Upon  termination  of  the  algorithm; 
will  serve  as  the  solution. 

During  each  epoch,  every  particle  is  accelerated 
towards  its  own  personal  best  as  well  as  in  the  direction 
of  the  global  best  position.  This  is  achieved  according 
to  the  following  two  expressions  [6]. 

vf=wXvf  +  c\Xrand]  -  Af) 

+  C2'><rand2()x{pf -Aj).  (8) 

Af  =  4  +  vf.  (9) 

where  k  -  U  •••♦  i  =  U  N,  c\  and  C2  are  two 
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positive  constants.  randl{  )  and  rand2{  )  are  two 
random  functions  in  the  range  [0,  1].  while  w  is  the 
inertia  weight.  Expression  (8)  is  used  to  calculate  the 
particle's  new  velocity  according  to  its  previous 
velocity  and  the  distances  of  its  current  position  from 
its  own  best  experience  (position)  and  the  group’s  best 
experience.  If  the  velocity  is  higher  than  a  certain  limit. 
Vntaxy  this  limit  will  be  used  as  the  new  velocity  for  this 
particle  in  this  dimension.  Thus,  keeping  the  particles 
within  the  search  space.  Then,  the  particle  flies  toward 
a  new  position  according  to  expression  (9). 

A  number  of  factors  affect  the  performance  of  the 
PSO.  First,  the  number  of  particles  in  the  swarm,  M, 
affects  the  runtime  significantly.  A  balance  between 
variety  (more  particles)  and  speed  (fewer  particles) 
must  be  considered.  Another  factor  is  the  maximum 
velocity  parameter,  v^ax  *  which  controls  the  maximum 
global  exploration  ability  PSO  can  have.  A  very  large 
value  for  this  parameter  can  result  in  oscillation.  On  the 
other  hand,  a  small  value  can  cause  the  particle  to 
become  trapped  in  local  minima.  The  inertia  weight,  w, 
is  employed  to  control  the  impact  of  the  previous 
history  of  velocities  on  the  current  velocity.  Thus,  it 
influences  the  trade-off  between  global  and  local 
exploration  abilities  of  the  particles.  A  larger  inertia 
weight  facilitates  global  exploration  and  the  search  of 
new  areas.  A  smaller  inertia  weight  tends  to  facilitate 
local  exploration  to  fine-tune  the  current  search  area. 
Suitable  selection  of  the  inertia  weight  can  provide  a 
balance  between  global  and  local  exploration  abilities 
[6]. 

For  the  implementation  used  in  this  paper,  the 
number  of  particles,  Af,  has  been  set  equal  to  35.  The 
maximum  velocity,  v^m  is  set  to  a  constant  value, 
which  is  equal  to  half  the  size  of  the  search  domain. 
The  inertia  weight  is  used  to  attenuate  the  magnitude  of 
the  velocity  updates  over  time.  This  attenuation  is  a 
linear  function  of  the  current  epoch  number.  Thus,  w 
linearly  decays  from  about  0.52  to  0.48,  The  constants 
Cl  and  C2  are  set  to  the  default  value  of  2. 


IV,  NUMERICAL  IMPLEMENTATION  AND 
SIMULATION  RESULTS 


The  proposed  approach  has  been  implemented  and 
computations  were  carried  out  for  three  different 
electromagnetic  device  configurations.  Throughout  the 
simulations.  comparisons  were  made  with 
computational  results  obtained  using  the  Quickfield  FE 
package  (Version  4.3)  and  the  core  nonlinear  B-H 
relation  was  assumed  as  shown  in  Fig.  1.  As  can  be 
seen  from  the  same  figure,  this  relation  was  reasonably 
approximated  using  expression  (3)  by  choosing  n  =  3 

and  C  =  i 


First,  simulations  were  carried  out,  using  the 
proposed  PSO  approach,  for  an  electromagnet  having 
length,  width,  depth  and  air-gap  length  of  l.Om,  1.0m, 
0.25m  and  0.1m,  respectively.  Magnetic  field  and 
current  density  distributions  were  investigated  subject 
to  coil  excitations  corresponding  to  current  densities  of 

1.5x10®  A/m^  and  0.75x10®  A/m^ .  These 
excitations  were  chosen  to  drive  the  core  in  the  initially 
linear  and  nonlinear  magnetization  ranges.  Beside 
giving  an  idea  about  the  relative  dimensions  of  the  coil 
in  comparison  with  the  core.  Figs.  2-5  demonsu^te  the 
simulation  results.  In  these  figures  a  vector  is  plotted  at 
the  center  of  every  triangular  sub-domain.  In  other 
words,  the  figures  give  also  information  about  the 
discretization  scheme. 

In  order  to  check  the  accuracy  of  the  proposed 
approach,  computations  were  compared  to  FE  identical 
simulations.  This  comparison  revealed  good  qualitative 
and  quantitative  agreement.  For  instance,  maximum 
flux  density  obtained  using  the  FE  approach  for  the 
high  and  low  excitation  levels  were  found  to  be 
1.13T  and  0.59  T,  respectively.  Furthermore,  Figs.  6- 
7  demonstrate  the  extent  to  which  results  obtained  using 
the  proposed  approach  match  FE  computations. 


Fig.l.  Assumed  exact  B^H  properties  and  its 
approximation  using  expression  (3). 


Fig.2,  Flux  density  distribution  obtained  using  the 
proposed  PSO  approach  for  the  electromagnet  subject 

to  an  excitation  of  1.5x10^  A/m^  (^max  =I  08T) . 
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Fig.3.  Magnetic  field  distribution  obtained  using  the 
proposed  PSO  approach  for  the  electromagnet  subject 
6  2 

to  an  excitation  of  1 .5x10  A/m  . 


(m) 


Fig.6.  Comparison  between  PSO  and  FE  computations 
for  the  particular  horizontal  contour  line  located  at  the 
center  of  the  electromagnet  window. 


(m) 


Fig.7.  Comparison  between  PSO  and  FE  computations 
for  the  particular  vertical  contour  line  located  at  the 
center  of  the  electromagnet  window. 


Fig.4.  Flux  density  distribution  obtained  using  the 
proposed  PSO  approach  for  the  electromagnet  subject 

to  an  excitation  of  0.75x10^  A/m^  (^max  =0.55T) . 


Fig.5.  Magnetic  field  distribution  obtained  using  the 
proposed  PSO  approach  for  the  electromagnet  subject 

to  an  excitation  of  0.75x10^  A/m^  . 


Second,  simulations  were  carried  out  for  an 
electromagnetic  actuator  having  window  dimensions, 
depth,  core  opening,  and  plunger  dimensions  of 
0.8mx0.9m,  0.2m,  0.50m  and  1.0mx0.4m,  respectively. 
Once  more,  magnetic  field  and  current  density 
distributions  were  investigated  for;  case#l  when  the 
plunger  is  in  contact  with  the  core  and,  case#2  when  it 
is  0.1m  apart.  Excitation  was  kept  constant 

5  2 

corresponding  to  a  cuaent  density  of  3x10  A/m  . 
Due  to  the  symmetry  of  this  configuration,  results  of 
these  distributions  are  only  given  in  the  right  half  of  the 
solution  domain  as  demonstrated  by  Figs.  8-11.  Again, 
these  figures  give  some  information  about  the  adopted 
discretization  scheme  as  well  as  the  relative  coil  to  core 
dimensions.  Figs.  12-13,  on  the  other  hand,  reveal  the 
good  agreement  between  results  obtained  using  the 
proposed  PSO  approach  and  those  obtained  through  the 
FE  technique.  It  should  also  be  mentioned  that,  using 
the  FE  technique,  the  maximum  flux  density  values 
were  found  to  be  1.16T  for  case#l  and  0,311  for 
case#2. 
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Fig.8.  Flux  density  distribution  obtained  using  the 
proposed  PSO  approach  for  the  electromagnetic 
actuator  when  the  plunger  is  in  contact  with  the  core 
(5n,a^=U2T). 


Fig.ll.  Magnetic  field  distribution  obtained  using  the 
proposed  PSO  approach  for  the  electromagnetic 
actuator  when  the  plunger  is  0.  Im  apart  from  the  core. 


Fig.9.  Magnetic  field  distribution  obtained  using  the 
proposed  PSO  approach  for  the  electromagnetic 
actuator  when  the  plunger  is  in  contact  with  the  core. 


Fig.  10.  Flux  density  distribution  obtained  using  the 
proposed  PSO  approach  for  the  electromagnetic 
actuator  when  the  plunger  is  0.1m  apart  from  the  core 
(5max  =  0.296  T). 


Fig.  12.  Comparison  between  PSO  and  FE  computations 
for  the  particular  horizontal  contour  line  located  at  the 
center  of  the  lower  core  yoke. 


Fig.  13.  Comparison  between  PSO  and  FE 
computations  for  the  particular  vertical  contour  line 
located  at  the  center  line  of  the  actuator. 
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Finally,  simulations  were  carried  out  for  a  switched 
reluctance  motor  (SRM)  whose  data  is  given  in  Table  I. 
Particular  attention  was  given  to  the  flux  density  and 
field  distributions  at  the  beginning  and  end  of  a 
stepping  phase.  Such  distributions,  driven  by  a  pair  of 
coil  excitation  corresponding  to  a  current  density  of 

3x10^  A/m^ ,  are  given  in  Figs.  14-17.  Obviously, 
these  figures  give  some  information  about  the  adopted 
discretization  scheme.  In  these  figures,  only  steady  state 
analysis  was  considered  and  no  motion  effects  were 
incorporated.  In  addition,  the  flux  was  assumed  to  be 
confined  inside  the  motor  (i.e.,  no  fiux  is  assumed  to 
cross  the  outer  stator  and  inner  rotor  diameters. 
Additional  comparisons  with  results  obtained  using  the 
FE  technique  are  given  in  Figs.  18  and  19.  It  is  also 
worth  mentioning  that  the  maximum  flux  density  values 
computed  using  the  FE  technique  at  the  beginning  and 
end  of  the  stepping  phase  were  found  to  be  0.08  T  and 
1.08  T,  rc.spcctivcly. 


TABLE I 

Data  of  the  simulated  Switched  Reluctance  Motor 


Number  of  Stator  Poles 

8 

Number  of  Stator  Poles 

6 

Stator  Pole  Arc 

15® 

Rotor  Pole  Arc 

15“ 

Stator  Core  Outer  Diameter 

0.100  m 

Stator  Core  Inner  Diameter 

0.750  m 

Diameter  at  Stator  Pole 

0.450  m 

Diameter  at  Rotor  Pole 

0.445  m 

Rotor  Core  Outer  Diameter 

0.250  m 

Rotor  Core  Inner  Diameter 

0.150  m 

Fig.  14.  Flux  density  distribution  obtained  using  the 
proposed  PSO  approach  for  the  SRM  at  the  beginning 
of  a  stepping  phase  ==  ^-07  T) . 


Fig.  15.  Magnetic  field  distribution  obtained  using  the 
proposed  PSO  approach  for  the  SRM  at  the  beginning 
of  a  stepping  phase. 


Fig.  16.  Flux  density  distribution  obtained  using  the 
proposed  PSO  approach  for  the  SRM  at  the  end  of  a 
stepping  phase  (B^gx  =1-13T). 


Fig.  17.  Magnetic  field  distribution  obtained  using  the 
proposed  PSO  approach  for  the  SRM  at  the  end  of  a 
stepping  phase. 
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Fig.  18.  Comparison  between  PSO  and  FE  computations 
for  the  contour  line  passing  through  the  centers  of  the 
two  energized  stator  poles  at  the  beginning  of  the 
stepping  phase. 


7.5  97.5  187.5 

Angle  (degrees) 

Fig.  19.  Comparison  between  PSO  and  FE  computations 
for  the  contour  arc  located  in  the  center  of  the  air-gap 
and  subtending  one  energized  stator  pole  to  the  other  at 
the  end  of  the  stepping  phase. 


IV.  DISCUSSION  AND  CONCLUSIONS 

In  light  of  the  proposed  PSO  approach  as  well  as  the 
experience  gained  while  performing  the  presented 
simulations,  the  following  few  remarks  should  be 
stressed: 

a)  Depending  on  the  problem  nature  and  number  of 
unknowns,  suitable  values  for  the  number  of 
particles  M  and  inertia  weight  w  should  be  first 
investigated.  In  our  case  it  took  some  time  to 
decide  upon  using  the  reported  M  and  w  values. 
Once  this  sort  of  tuning  is  achieved,  the  proposed 
approach  gives  reasonable  results  in  computational 
time  comparable  to  that  of  the  FE  technique. 

b)  It  can  be  observed  from  the  presented  results  that 
good  agreement  with  FE  computations  can  be 
achieved  using  the  proposed  approach.  It  should  be 
pointed  out,  however,  that  some  of  the 
discrepancies  between  results  of  both  techniques 


stem  from  the  non-identical  discretization  schemes, 

c)  An  important  feature  of  the  proposed  search 
approach  is  its  simplicity  and  ability  to  handle 
adjacent  discretizations  of  dramatically  different 
dimensions  without  the  risk  of  running  into 
numerical  difficulties  related  to  matrix  inversion. 
This  further  highlights  the  computational  memory- 
wise  efficiency  of  the  approach. 

More  efforts  are  planned  to  further  investigate  the 
proposed  approach  in  different  time  harmonic 
problems.  It  should  also  be  stressed  that  PSO  has  been 
previously  used  in  device  dimension  optimization  (see, 
for  example,  [7]).  Coupling  this  capability  to  the 
proposed  field  analysis  approach  might  pave  the  way 
towards  the  ability  to  simultaneously  optimize  the 
dimensions  and  compute  fields  in  devices  involving 
nonlinear  media. 


REFERENCES 

[1]  J.K.  Sykulski,  Computational  Magnetics^  Chapman 
&  Hall,  London,  1995. 

[2]  A.A.  Adly  and  S.K.  Abd-El-Hafiz,  "Automated 
two-dimensional  field  computation  in  nonlinear 
magnetic  media  using  Hopfield  neural  networks,” 
IEEE  Trans.  Magn.,  Vol.  38,  pp.  2364-2366, 2002. 

[3]  M.V.K.  Chari  and  S.J.  Salon.  Numerical  Methods 
in  Electromagnetism,  Academic  Press,  San  Diego 
CA,  2000. 

[4]  A.A.  Adly  and  S.K.  Abd-El-Hafiz,  "Utilizing 
Hopfield  neural  networks  in  analysis  of  reluctance 
motors,"  IEEE  Trans.  Magn.,  vol.  36,  pp.  3147- 
3149, 2000. 

[5]  J.  Kennedy  and  R.  Eberhart,  "Particle  Swarm 
Optimization,"  Proceedings  of  the  IEEE 
International  Conference  on  Neural  Networks, 
Perth,  Australia,  pp.  1942-1948, 1995. 

[6]  Y.  Shi  and  R.  Eberhart,  "A  Modifred  Particle 
Swarm  Optimizer,”  Proceedings  of  the  IEEE 
International  Conference  on  Evolutionary 
Computation,  Anchorage,  Alaska,  pp,  69-73, 1998. 

[7]  G.  Ciuprina,  D.  loan  and  L'Munteanu,  "Use  of 
Intelligent-Particle  Swarm  Optimization  in 
Electromagnetics,"  IEEE  Trans.  Magn.,  vol.  38, 
pp.  1037-1040.2002. 


Amr  A.  Adly  received  the  BS 
and  MS  degrees  in  Electrical 
Engineering  from  Cairo 
University,  Egypt,  in  1984  and 
1987,  respectively,  and  the  PhD 
degree  in  Electrical  Engineering 
from  the  University  of 
Maryland,  College  Park, 
Maryland,  USA,  in  1992.  In 


209 


ACES  JOURNAL,  VOL.  18.  NO.  3,  NOVEMBER  2003 


1993  he  joined  LDJ  Electronics,  Troy,  Michigan,  USA, 
as  a  Senior  Engineer/Scientist.  In  1994,  he  was 
appointed  as  an  Assistant  Professor  at  the  Electrical 
Power  and  Machines  Department,  Faculty  of 
Engineering,  Cairo  University,  Giza,  Egypt.  Since 
1999,  he  has  been  working  as  an  Associate  Professor  in 
the  same  department.  His  research  interests  include 
electromagnetics,  modeling  and  simulation  of  complex 
magnetic  materials,  electrical  machines, 
superconductivity,  magnetic  recording,  magnetic 
measurement  instrumentation,  and  magneto-hydro¬ 
dynamics.  In  1994,  he  received  the  Egyptian  State  Prize 
for  achievements  in  Engineering  Sciences.  He  has 
authored  more  than  60  journal  and  refereed  conference 
papers.  Dr.  Adly  is  a  senior  member  of  the  IEEE 
Magnetics  Society  and  a  member  of  the  ACES  Society. 

SSalwa  K.  Abd-El-Hafiz 
received  the  BS  degree  in 
Electrical  Engineering  from 
Cairo  University,  Egypt,  in  1986 
and  the  MS  and  PhD  degrees  in 
Computer  Science  from  the 
University  of  Maryland,  College 
Park,  Maryland,  USA,  in  1990 
and  1994,  respectively.  In  1994,  she  was  appointed  as 
an  Assistant  Professor  at  the  Engineering  Mathematics 
Department,  Faculty  of  Engineering,  Cairo  University, 
Giza,  Egypt.  Since  1999,  she  has  been  working  as  an 
Associate  Professor  in  the  same  department.  Her 
research  interests  include  software  analysis  and 
specification,  software  measurements,  software  reverse 
engineering,  neural  networks,  evolutionary  computation 
techniques,  knowledge-based  systems  and  numerical 
analysis.  In  2001,  she  received  the  Egyptian  State  Prize 
for  achievements  in  Engineering  Sciences.  Dr.  Abd-El- 
Hafiz  is  a  member  of  the  IEEE  Software  Computer 
Society. 


ACES  JOURNAL,  VOL.  18,  NO.  3,  NOVEMBER  2003 


AN  APPROACH  TO  MULTI-RESOLUTION  IN  TIME  DOMAIN  BASED  ON  THE  DISCRETE 

WAVELET  TRANSFORM 


C.  Represa(*),  C.  Pereira(*),  A.C.L.  Cabeceira  (♦*),  I.  Barba  (♦*),  J.  Represa  (**) 

(*)  Dpt.  of  Electromechanical  Engineering.  University  of  Burgos.  09001  Burgos.  Spain. 

Phone:  +34  947258830  FAX:  +34  94725883 1  E-Mail:  creDresa@ubu.es 
(**)  Dpt.  of  Electricity  and  Electron.  University  of  Valladolid.  4701 1  Valladolid.  Spain. 
Phone:  +34  983423223  FAX:  +34  983423225  E-Mail:  ibarba@ee.uva.es 


Abstract^-  In  this  paper,  an  approach  to  multi- 
resolution  in  time  domain  (MRTD)  is  presented. 
Maxwell  equations  are  discretized  using  finite 
differences  in  time  and  a  derivative  matrix  in 
space  that  allows  any  desired  level  of  spatial 
resolution.  This  derivative  matrix  acts  on  the 
coefficients  that  represent  the  expansion  of  the 
^eld  components.  These  coefficients  are 
calculated  by  means  of  the  Discrete  Wavelet 
Transforms  (DWT).  In  this  work  hard  (PEC 
and  PMC)  boundary  conditions  have  been 
introduced  into  the  algorithm  using  the  method 
of  images.  This  approach  Is  valid  for  any  kind  of 
wavelet  functions.  Stability  and  dispersion 
properties  are  also  investigated.  Some  numerical 
results,  showing  multi-resolution  properties  are 
presented. 

Keywords.-  Multi-resolution  in  time  domain, 
MRTD,  Daubechies  wavelets,  Discrete  Wavelet 
Transform,  DWT. 


1.-  Introduction 

The  development  of  electromagnetic  fields  in  scale 
and  wavelet  functions  [1],  has  given  rise  to  the 
techniques  known  as  Multi  Resolution  in  Time 
Domain  (MRTD).  These  techniques  are  based  on 
the  possibility  of  increasing  the  resolution  of  a 
signal,  from  a  coarse  level  to  a  fine  one,  using  low- 
resolution  functions  (scale  functions),  combined 
with  others  of  intermediate  resolution  levels 
(wavelet  functions).  Different  functions  have  been 
used  in  this  type  of  analysis:  Battle-Lemarie  [2], 
Haar  [3]  or  Daubechies  [4]. 

In  these  techniques,  the  electric  and  magnetic  field 
components  of  Maxwell's  curl  equations  are 
represented  by  a  manyfold  expansion  in  scale  and 
wavelet  functions  with  respect  to  space,  and  step 
functions  with  respect  to  time.  The  method  of 


‘  This  work  was  sponsored  in  part  by  the  Spanish  Ministry  of 
Science  and  Technology,  under  its  project  TIC-2000- 1612-C03- 
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moments  [5]  allows  to  obtain  a  set  of  equations 
similar  to  the  one  used  in  FDTD  (identical,  in  the 
case  of  Haar  functions,  to  the  Yee’s  algorithm). 
Instead  of  the  field  values,  the  coefficients  of  this 
expansion  are  used  in  such  scheme. 

One  of  the  points  that  affect  the  complexity  of 
these  schemes  is  the  type  of  functions  used  for  the 
expansion.  For  example,  if  Battle-Lemari6 
functions  are  used,  an  infinite  number  of 
coefficients  must  be  used,  because  of  the  non¬ 
compact  support  of  these  fimctions.  Actually,  we 
can  truncate  the  expansion  and  use  a  reduced 
number  of  them,  obtaining  an  approximate 
solution.  This  problem  does  not  appear  when,  for 
example,  Haar  functions  are  used,  since  they  are 
compactly  supported,  so  it  is  not  necessary  to 
trancate  the  expansion. 

A  second  point  is  the  desired  resolution  level  for 
the  solution  of  the  problem.  The  above  mentioned 
approaches  [2-4]  use  scale  functions  for  a  first 
approximation  of  the  fields;  in  this  case  simple 
functions  are  obtained,  but,  if  an  increasing  of  the 
resolution  is  desired,  adding  different  levels  of 
wavelet  functions,  the  solution  becomes  unfeasible, 
as  the  scheme  must  be  modified  for  every  new 
level  we  add.  This  problem  can  be  solved  if  a 
Discrete  Wavelet  Transform  (DWT)  is  used  to 
derive  the  coefficients  of  the  expansions,  as  it  h^ 
already  been  used  4o  obtain  the  time  domain 
solution  of  electrical  networks  [6]  - 

In  this  work,  compactly  supported  Daubechies 
wavelet  functions  [1]  have  been  used  to  make  an 
expansion  of  the  fields.  Its  coefficients  are 
computed  through  a  Discrete  Wavelet  Transform. 
A  derivative  matrix,  also  obtained  by  means  of  a 
DWT,  allows  us  to  compute  the  value  of  the 
electric  and  magnetic  fields  at  desired  spatial 
positions  and  with  the  desired  resolution.  We  have 
found  the  stability  criterion  of  the  algorithm  for 
different  wavelet  functions  and  in  some  cases  it  is 
wider  than  the  FDTD  stability  criterion  for  the 
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same  spatial  resolution  (twice  in  the  case  of  Haar 
wavelet  functions). 


2.-  Formulation 


The  simplest  case  of  a  TEM  plane  wave 
propagating  in  a  homogeneous,  linear,  isotropic  and 
non-dispersive  media,  with  fields  Ex  and  Hy  is 
analyzed.  Then,  the  equations  to  solve  are: 


dH, 

^  dt 

y 

az  “  ^  at 


dz 

aH, 


(1) 


In  the  proposed  scheme,  every  field  component  is 
expanded  in  terms  of  scaling  and  wavelet  fonctions, 
as  shown  in  equation  (2),  where  the  coefficients  of 
the  lower  resolution  level  are  given  by  the 
scaling  function  ((j))  and  the  coefficients  of  the 
successive  resolution  level  (6/)  are  given  by  the 
wavelet  functions  (\|/): 


I"a:K(z) 

" «  (2) 

+  2:i"bW(z) 

,  k  i=0 

The  spatial  discretization  is  achieved  by  dividing 
the  domain  into  N  cells  of  size  A,  as  a  starting 
point,  then  getting  N  sampling  points  in  the  center 
of  each  cell.  This  is  what  we  call  resolution  level  0 
(j  =  0).  If  each  cell  is  now  divided  into  2^  points, 
the  maximum  level  of  desired  resolution  J  (j  =  J)  is 
reached  and  then  the  simulation  domain  has 


F''(z)  =  X''ai<t.;!(z)  = 

k 


sampling  points  with  a  distance  of  Az  between 
them.  This  is  illustrated  in  figure  1 . 

I  ■  I  •  j  J-o 


1  •  1  •  1  •  1  •_ 

•  l^_l 

i- 

- -  A  - 

LazH 

Fig.  L-  Spatial  discretization  for  J-2. 


The  coefficients  are  enough  to  represent  the  real 
values  of  the  field  at  each  sampling  point  and  they 
are  the  initial  and  the  final  step  of  the 
decomposition  process  involved  in  the  Discrete 
Wavelet  Transform.  The  multi-resolution  analysis 
performed  with  the  wavelet  transform  can  be 
understood  as  a  digital  filtering  process  where  a 
signal  is  decomposed  in  two  parts,  one  containing 
the  low  frequencies,  the  scale  functions  (f),  and 
other  part  containing  the  high  frequencies,  the 
wavelet  functions  (^.  The  multi-resolution  - 
representation  through  the  Discrete  Wavelet 
Transform  (DWT)  is  provided  by  successive  filter 
banks  stages,  each  one  containing  a  low-pass  and  a 
high-pass  filter,  described  in  terms  of  the 
coefficients  of  their  impulse  responses  L(m)  and 
H(m)  (meZ)  respectively  [1].  The  filtering  process 
means  successive  convolutions  of  the  field  and  the 
filter  coefficients  followed  by  a  decimation  process 
that  retains  only  the  even  indexes.  The  inverse 
transform  (IDWT)  consists  in  successive 
interpolations  followed  by  convolutions  that  gives 
the  field  values  using  the  coefficients  of  the 
previous  decomposition.  This  is  sketched  in  figure 
2. 


D^composilion  Wave/ef  CoefUdents  Reconsinidhn 

EWT  IDWT 

Fig.  2.-  A  step  of  the  filtering  process:  the  low  and  high-pass  filters  outputs  the  coefficients  d  and  b! 
respectively. 

Then,  eqs.  (1)  arc  discretized  using  centered  finite  the  other  field  coefficients  at  an  intermediate  time 

differences  for  the  time  derivatives  and  a  spatial  step: 

derivative  operator  computed  through  a  DWT 

matrix  for  the  spatial  derivatives.  Each  coefficient 

is,  then,  obtained  from  its  values  at  the  same  point 

in  a  previous  time  step  and  from  the  derivation  of 
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-DVe)" 

glHl  _  p 


(3) 


where  At  and  A  are  the  time  and  space  steps 
respectively,  and  D**  represents  the  spatial 
derivative  operator  of  level  J.  This  operator  [7]  acts 
on  the  coefficients  a  and  b  of  the  field  transform  (3) 
and  it  is  denoted  by  a  matrix  D  [8]  which  can  be 
computed  previously.  In  this  way,  eqs.  (3)  can  be 
rewritten  as  follows  for  an  arbitrary  level  of 
resolution  J: 


J  n 

r  |-|J  ft+l 

E’'  Atl  ,  ’ 

"  J 

E**. 

E**  e  A 

H*'_ 

[H*  J  M  A  1 

(4) 


where  This  matrix  D  is  banded  with  a 

limited  band  width,  due  to  the  compact  support  of 
the  Daubechies  wavelets,  and  this  width  increases 
as  the  size  of  the  filters  does.  Its  elements  are 
integrals  of  the  scaling  and  wavelet  functions  and 
its  discrete  spatial  translations  of  the  first  derivative 
(eqs.  (6.1)-(6.4)).  Because  of  multi-resolution,  this 
matrix  can  be  decomposed  in  a  set  of  submatrices 
of  lower  resolution,  so  if  we  have  named  the 
original  matrix  of  highest  resolution  O'*,  it  is 
constructed  with  the  one  level  lower  matrices  A,  B, 
F,  D,  then: 


(r*-' 

where A<  =  ^i}.BU|j3j}.ri4i},D‘ 41}. 
The  components  of  every  matrix  can  be  calculated 
in  a  recursive  way,  so: 

a||=2'ai^  (6) 

and  similar  expressions  for  the  other  coefficients,  p, 
y  and  d.  The  expression  aj  is  obtained  by  means  of 
the  following  expression: 

a,  =  J\|/(z-l)~»)/(z)dz  (7) 

Similar  expressions,  combining  scale  and  wavelet 
functions  can  be  used  to  calculate  the  other 
coefficients  p  (scale, wavelet),  y  (wavelet  and  scale) 
and  d  (scale  and  scale).  The  scale  and  wavelet 
functions  in  every  case  are  obtained  by  means  of 
the  low-pass  and  the  high  pass  filters,  L(m)  and 
H(m)  respectively: 


(t)(t)  =  2L(m)-y24.(2t-m)  (8) 

in=0 

v{t)  =  |;H(m),^<K2t-m)  (9) 

m=o 

being  Z/the  filter  length.  In  figure  3  three  different 
aspects  of  the  derivative  matrix  D  for  three  levels 
of  resolution  are  depicted,  where  only  the  nonzero 
elements  have  been  plotted. 


(c) 

'  Fig.  i.-  piree  different  aspects  of  the  derivative 
matrix  if  for  three  resolution  levels:  (a)  j-0,  (b) 
7=/,  (c)j^2. 


The  use  of  this  derivative  matrix  allows  us  to 
achieve  any  level  of  resolution  without 
modifications  of  the  algorithm.  Moreover,  we  can 
choose  at  the  beginning  of  the  simulation  what  type 
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of  wavelet  function  is  going  to  be  used,  and  what 
regions  of  the  simulation  domain  we  want  tosolve 
at  higher  resolution.  This  feature  will  be  shown 
further  in  the  examples. 


3.-  Boundary  conditions. 

The  derivative  matrix  D  is  built-up  in  a  cyclic  form 
in  such  a  way  that  it  treats  the  coefficients  on  both 
boundaries  as  contiguous.  It  results  in  a  cyclic 
space,  that  is,  an  infinite  unbounded  space.  As  the 
simulation  domain  is  not  cyclic  at  all,  the  character 
of  the  matrix  must  be  modified.  To  do  so,  we 
propose  to  add  adjacent  columns  with  the  elements 
related  to  the  boundaries,  and  getting  an  extended 
matrix;  i.e. 


'a,  a.,  0  a,' 

a,  Oo  a_,  0 

0  a,  a,  a., 
^a_,  0  a, 


KL 


Uo 

a-i 

0 

0 

0  ' 

0 

i«. 

a., 

0 

0 

0 

1  0 

a, 

Cto 

a-i 

1  0 

0 

1  0 

0 

a. 

oto 

|a.i> 

and  similar  expressions  for  },{/},  and  {d}.  The 
matrix  D  acting  on  the  column  vector  (b'*,  a"*)  at 
resolution  level  J  implies  that  its  extension  must  act 
on  an  extended  vector  too.  The  additional 
coefficients  needed  to  represent  the  boundaries  can 
be  obtained  using  the  method  of  images  for  PEC  or 
PMC  walls.  If  di  and  h/i  arc  the  scale  and  wavelet 
coefficients  respectively,  belonging  to  a  border  cell 
collocated  at  position  we  obtain  the 

additional  coefficients  a*  and  b'  this  way: 


More  complex  structures  require  more 
sophisticated  techniques  not  studied  here. 

4.-  Stability  analysis 

For  the  algorithm  to  be  stable,  the  pair  At  and  Az 
must  be  chosen  in  a  correct  way.  We  choose  these, 
values  following  the  derivation  given  in  [9]  where 
the  stability  problem  is  treated  as  an  eigenvalue 
problem.  This  means  that  plane  wave  cigenmodes 
will  be  assumed  to  propagate  in  the  numerical  data 
space.  The  spectrum  of  eigenvalues  for  these 
modes  due  to  the  numerical  space  differentiation 
process  will  be  determined  and  compared  to  the 
stable  spectrum  of  eigenvalues  determined  by  the 
numerical  time  differentiation  process.  By 
requiring  the  complete  spectrum  of  spatial 
eigenvalues  to  be  contained  within  the  stable  range, 
it  is  ensured  that  all  possible  numerical  wave 
modes  in  the  grid  are  stable.  In  this  way,  equations 
(I)  can  then  be  rewritten  as  follow: 


At 

E''"*  —  E*' 
At 


(10) 


and  this  equations  can  be  split  into  two  eigenvalue 
problems  concerning  time  and  space: 


H"*' 

At 

^n+1  n 

At 


=  AH" 


=  AE" 


(H) 


1 

p-Az 


I 


1 

e- Az 


£d,H:;UAE,'^ 


(12) 


•  for  even  symmetry: 


•  for  odd  symmetry: 


on  the  ri^ 


on  the  left 


on  the  rigth 


on  the  left 


These  boundary  conditions  allow  us  to  implement 
in  a  easy  way  infinite  electric  or  magnetic  walls. 


where  only  .scaling  functions  of  resolution  level 
J=0  have  been  used.  To  avofd  having  any  mode 
--  “increasing  without  limit  during  normal  time¬ 
stepping  it  is  found  from  equation  (11)  that  the 
stable  spectrum  of  eigenvalues  is: 

|lm(A)|s2/At  (13) 

In  order  to  solve  equations  (12)  we  introduce  a 
typical  mode  of  the  spatial  spectrum  like  (14)  into 
the  equations: 

"E  =  E 

'  °  ^  (14) 

Then  we  get  the  following  set  of  equations: 
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H-AZ  Y  / 


ai(««l-A2)ni(-X+A2)  . 


Simplifying  and  using  Euler's  identity  we  get  that 
the  spatial  eigenvalues  are  given  by: 


lm(A)  =  —  y  d,sen(~klA2)  (16) 

AzT 

Therefore,  the  maximum  value  of  this  eigenvalues 


|lm(A)|<i£|d,|  (17) 

To  guarantee  numerical  stability,  the  range  of 
spatial  modes  must  be  contained  completely  within 
the  stable  range  of  time^stepping  eigenvalues  set  by 
(17)  and  so: 

<“> 


Therefore,  the  upper  bound  of  the  time  step  AtwajD 


Defining  a  stability  factor  like  (20): 


propagation,  that  is,  the  numerical  phase  velocity 
given  by  the  algoridim  differ  from  the  phase 
velocity  of  the  wave  in  the  medium  [11].  This 
fictitious  dispersive  behaviour  must  be  taken  into 
account  specially  in  considering  large  structures  of 
simulation  because  significant  differences  between 
real  and  numerical  phase  can  be  obtained.  Some 
studies  about  dispersion  properties  of  MRTD 
methods  based  on  Galerkin  procedures  [12],  [13] 
have  been  done.  Now  we  proceed  to  study  the 
dispersion  characteristics  of  the  MRTD  algorithm 
based  on  the  DWT  using  Daubechies'  wavelets 
functions.  To  do  that,  a  plane  monochromatic 
travelling-wave  (22)  is  introduced  into  the 
discretized  Maxwell  equations  (4)  and  then  we 
search  for  the  relationship  between  the  angular 
frequency  (o  and  the  numerical  wave  number  k : 

npC  -.-C  ftjteAz-wriAt) 

As  starting  point,  we  expand  the  field  components 
using  only  scaling  functions  of  resolution  level 
J=0,  so  ^  =  (1)  and  the  derivative  matrix  D®  will  be 
composed  of  elements  The  discretized 
Maxwell  equations  can  be  written  in  this  way: 

E  Z1Z[_  I  _ 

Substituting  (22)  into  (23),  simplifying  and  using 
Euler's  identity  we  get  this  set  of  equations: 


equation  (19)  establishes  that  the  stability  factor 
must  be  contained  within  the  range : 

0<Si^  (21) 


“Depending  on  the  iype  of  wavelet  function  used, 
this  range  implies  that  the  time  step  can  be  set  to  a 
different  value  for  the  same  spatial  resolution.  If 
Haar  wavelet  functions  are  used,  this  value  can  be 
set  double  than  the  time  step  set  in  FDTD  for  the 
same  spatial  resolution. 


5.-  Dispersion  properties 

The  use  of  numerical  techniques  for  solving 
electromagnetic  problems  as  TLM  [10],  FDTD  or 
MRTD  require  always  a  discretization  process. 
Such  a  process  result  in  a  phase  error  of  the  field 


(oAt  At 


2EoSen*^  =  Xd,HoSen(-  klAzj 

y  \(24) 

(oAt  At  1  (  r..  1 

2HoSen— = — —  Y  djEoSeni-klAzI 
2,  p  AZl  I 

and  from  them  we  obtain  the  dispersion 
relationship  where  ds  the  speed  of  light  in  our 
medium: 

^sen-5^  =  2;d,sen(-klAz)  (25) 


This  expression  may  be  reduced  to  the  ordinary 
expression  co  =  ck ,  that  is,  the  non-dispersive  case, 
if  spatial  resolution  Az  is  very  small  in  comparison 
with  the  wavelength.  If  Haar  wavelet  functions  are 
used  then  we  also  obtain  a  non  dispersive  case  at 
the  stability  factor  s  =  2  (i.e.  the  maximum  stable 
value).  In  any  case,  dealing  with  broadband  signals, 
different  expressions  will  be  obtained  depending  on 
the  relation  between  the  wavelength  A  and  the 
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spatial  resolution  Az,  and  on  the  type  of  wavelet 
function  we  use.  In  figure  4  it  is  plotted  the 
nonmalized  phase  velocity  versus  spatial  resolution 
given  by  the  expresion  (25)  for  two  different  values 
of  the  stability  factor.  We  can  appreciate  the  feature 
mentioned  before,  that  is,  the  normalized  phase 
velocity  tend  to  unity  as  we  incresc  the  resolution. 
It  can  also  be  seen  that  the  higher  the  order  of  the 
wavelet  function,  the  better  dispersion 
characteristics. 


Normetaed  phsso  vetodty  for  S  =05 


Normaiizod  pheso  volodty  for  S  ■=  1 


(b) 

Fig.  4.-  Normalized  numerical  phase  velocity 
versus  spatial  resolution  in  lambda  units  for 
different  wavelet  functions: 

(a)  with  a  stability  factor  s-0,5,  and 

(b)  with  a  stability  factor  s= 1.0 


In  figure  5  it  is  depicted  a  gaussian  pulse 
propagation  using  different  types  of  wavelet 
functions:  (a)  Haar  wavelet  functions,  and  (b) 
Daubcchics  Di  wavelet  functions.  The  case  (b) 
exhibit  a  better  dispersive  behaviour  than  the  case 
(a)  as  it  can  be  expected  from  previous  figure  4. 


(b) 


Fig.  5.-  Three  snapshots  of  a  Gaussian  pulse 
propagation  using  (a)  Haar  wavelet  functions,  and 
(b)  Daubechies  Db2  functions,  showing  it 
dispersion  properties. 


6*-  Results 

To  validate  our  proposed  technique  some  aditional 
examples  have  been  performed.  First  of  all  we  have 
evaluated  the  multi-resolution  property  of  the 
algorithm.  To  do  that,  we  have  simulated  a 
gaussian  pulse  excitation  into  a  one-dimensional 
cavity  of  1  m  long.  The  simulation  domain  is  split 
into  two  parts  of  different  resolution,  one  of  high 
resolution  on  the  left  side  and  400  mm  long,  and 
other  part  of  low  resolution  on  the  right  side  and 
600  mm  long.  Figure  6  displays  this  simulation 
where  it  can  be  clearly  seen  the  differences 
between  the  two  parts.  Figure  7  shows  in  detail  the 
transition  between  zones. 
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Fig.  6.-  Gaussian  pulse  propagation  using  Haar 
wavelet  Junctions.  The  simulation  domain  is  split 
into  two  zones  of  different  resolution. 


Time  =  1450  ps 


Fig.  7.-  Details  of  the  transition  between  two  zones 
of  different  resolution. 


We  have  also  evaluated  the  resonant  frequencies  of 
a  one-dimensional  rectangular  cavity  with  a  100 
mm  distance  between  electric  walls  (PEC).  A 
spatial  resolution  of  Az  ~  0.625  mm  and  a  time 
discretization  of  At  =  3.127  ps  have  been  chosen 
(that  means,  that,  according  to  the  stability  criterion 
we  choose  s  =  1.5).  The  electric  and  magnetic 
fields  are  specified  at  ^0  and  at  /=l/2At 
respectively, 


where  the  width  of  the  pulse  is  =  10  and  its 
center  kc  =  80  in  terms  of  the  index  of  the  spatial 
mesh.  After  4096  time  steps  with  a  spectral 
resolution  of  Af  =  78.1  MHz,  the  resonance 
spectrum  obtained  have  been  plotted  in  figure  8 


Fig.  8.-  Resonant  frequencies  obtained  in  a  cavity 
of  length  100  mm  using  Haar  wavelet  Junctions. 
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7.-  Conclusions 

A  new  approach  to  Multi-Resolution  in  Time- 
Domain  has  been  investigated.  In  this  study  a 
derivative  matrix  is  used  to  calculate  the  spatial 
derivatives  of  the  electromagnetic  fields.  This 
matrix  acts  on  the  coefficients  of  the  wavelet 
expansions  of  the  fields  obtained  from  the  Discrete 
Wavelet  Transform.  The  use  of  this  matrix  allows 
us  to  solve  the  discretized  Maxwell  equations  at 
any  level  of  desired  spatial  resolution  and  wherever 
we  may  want.  It  has  been  shown  that  different 
types  of  Daubechies  wavelet  functions  exhibit 
better  dispersion  properties  for  the  same  spatial 
resolution  and  that  the  time  step  can  be  chosen 
bigger  than  in  other  time-domain  methods  with  the 
same  mesh  size.  In  this  way,  wc  can  also  choose  at 
the  beginning  of  the  simulation  different  types  of 
wavelet  functions  in  order  to  improve  our  results. 
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or  they  may  focus  on  specific  applications,  techniqcs, 
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physical  insight  or  understanding  is  achieved.  - 

4  New  computational  techniqcs  ,  or  new  applications  for 

existing  conputational  techniques  orcodes. 

5  Tricks  of  the  trade”  in  selecting  and  applying  codes 

and  techniques. 

6  New  codes, algorithms,  code  enhancement,  and  code 
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improvement.  Note:  Code  algorithm)  capability 
descriptions  are  not  acceptable,  unless  they  contain 
sufficient  technical  material  to  justi^^  consideration. 
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computer-aided  design)  or  in  output  (whether  it  be 
tabular,  graphical,  statistical,  Fourier-transformed,  or 
otherwise  signal-processed).  Material  dealing  with 
input/output  database  management,  output  interpretation, 
or  other  input/output  issues  will  also  be  considered  for 
publication. 

8  Computer  hardware  issues.  This  is  the  category  for 

analysis  of  hardware  capabilities  and  limitations  of 
various  types  of  electromagnetics  computational 
requirements.  Vector  and  parallel  computational 
techniques  and  implementation  are  of  particular  interest. 
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biological  effects,  electromagnetic  pulse  (EMP), 
electromagnetic  interference  (EMI),  electromagnetic 
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dielectric,  magnetic  and  nonlinear  materials,  microwave 
components,  MEMS  technology,  MMIC  technology,  remote 
sensing  and  geometrical  and  physical  optics,  radar  and 
communications  systems,  fiber  optics,  plasmas,  particle 
acceleratois,  generators  and  motors,  electromagnetic  wave 
propagation,  non-destructive  evaluation,  eddy  currents,  and 
inverse  scatterir^. 

Techniques  of  interest  include  frequency-domain  and  time- 
domain  techniques,  integral  equation  and  difErential  equation 
techniques,  diffraction  theories,  physical  optics,  moment 
methods,  finite  differences  and  finite  element  techniques, 
modal  expansions,  perturbation  methods,  and  hybrid  methods. 
This  list  is  not  exhaustive. 
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unsuccessful  efforts  in  applied  computational 
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means  to  discuss  problem  areas  in  electromagnetic  modeling. 
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results  in  computational  electromgnetics  will  be  considered 
for  publication  only  if  a  reasonable  expectation  of  success 
(and  a  reasonable  effort)  are  reflected.  Moreover,  such 
material  must  represent  a  problem  area  of  potential  interest  to 
the  ACES  membership. 

Where  possible  and  appropriate,  authors  are  required  to 
provide  statements  of  quantitative  accuracy  for  measured 
and/or  computed  data.  This  issue  is  discussed  in  “Accuracy 
&  Publication:  Requiring,  quantitative  accuracy  statements  to 
accompany  data,”  by  E.  K.  Miller,  ACES  Newsletter,  Vol.  9, 
No.  3,  pp.  23-29, 1994,  ISBN  1056^9170. 
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prepared  on  srandard  8.5cl  1  inch  paper. 
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